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ABSTRACT 


Numerous computational methods for generating and 
simulating binary and grey-level natural dig ital-image 
textures are proposed using a variety of stochastic 
models. Pictorial results of each method are given and 
various aspects of each approach are discussed. The 
quality of the natural texture simulations depends on the 
computation time for data collection, computation time for 
generation, and storage used in each process. In most 
cases, as computation time and data storage increase, the 
visual match between the texture simulation and the parent 
texture improves. Many textures are adequately simulated 
using simple models thus providing a potentially great 
information compression for many applications. 

In most of the texture synthesis methods presented in 
this thesis, pixel values are generated one-at-a-1 ime 
according to both the given model and the values of pixels 
previously generated in the synthesis until the image 
space is completely filled. Nth-order joint density 
functions estimated from a natural texture sample were 
used for this purpose in one method. The results are 
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excellent but the storage required, even for binary 
textures, is large. Therefore, a much simpler first-order 
linear, autoregressive model was applied to the texture 
synthesis problem. Using this model on both binary and 
continuous-tone textures, each pixel is generated as a 
linear combination of previously generated pixels plus 
stationary noise. The results indicate that many textures 
are satisfactorily simulated using this approach. 

By adding cross-product terms, the first-order linear 
model is extended to a second-order linear model. The 
simulation results improve slightly but the number of 
computations required for the statistics collection 
process increases drastically. Non-stationary noise was 
then used in the synthesis process and further 
improvements in the quality of the simulations are 
achieved at the cost of increased storage. 

Methods of texture simulation using more than one 
model are studied in this thesis. These multiple-model 
are useful for many textures, especially those with 
macro-structure. They also improve the fit of the model 
when applied to the parent texture data and therefore may 
produce improved simulations. 

A final model, called the best-fit model, generates 
texture simulations directly from the parent texture 
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itsel Each pixel in the synthesis image is generated 
based on the similarity of its previously-generated. 


neighboring pixel values to pixel values in all 
similarly-shaped neighborhoods in the parent texture. The 
measures of similarity at all points in the parent 
texture, along with a random variable, are used to 
generate the next pixel value in the synthesized image. 
The synthesis results using model are excellent but the 
synthesis process is very computationally demanding. 

Although the success of texture synthesis is highly 
dependent on the texture itself and the modeling method 
chosen, general conclusions regarding the performance of 
various techniques are given. Methods of texture 
segmentation and identification based on texture synthesis 
results are also presented. 
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CHAPTER 1 


INTRODUCTION 


1.1 Introduction 

Texture is important characteristic for the analysis 
of many types of images. It is an important feature for 
discrimination and identification of regions in images and 
as a result the vast majority of work on texture has 
concentrated on these applications. Although many 
different texture discrimination techniques have been 
developed, most are ad hoc . 

The problem inverse to texture analysis is texture 
synthesis, or the generation of image fields having 
analytical and visual characteristics similar to natural 
textures. Texture synthesis has been over-shadowed by the 
emphasis placed on the discrimination problem and its 
applications. Little work has been done on synthesis even 
though numerous applications exist. For example, 
intelligent image sensors could transmit boundaries of 
textured image regions. Based on statistics gathered by 
the sensor, this region could be reconstructed using 
simulation techniques with little or no loss of 
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information. The result is excellent data compression. 


Texture synthesis can also be used as a texture 
analysis tool leading to a better understanding of 
textures and their perception by humans as v/ell as 
improved methods of discrimination. By carefully 
controlling the statistics of a texture in a synthesis 
process visual changes in texture are observed. Thus, 
texture synthesis methods allow researchers to identify 
and measure the information content of individual 
statistical measurements. By assembling these 
measurements and incorporating them into a texture 
simulation process, statistics may be measured from a 
parent texture and used to produce a texture simulation. 
The degree to which the parent and simulation are visually 
similar indicates the value of the statistical 
measurements and the model used in the simulation process. 
Given a group of statistical measurements which are 
proposed to be useful texture measures, possibly the best 
may be chosen based on the quality of the corresponding 
texture simulations. In this way, researchers are able to 
develop better discrim ination as well as better synthesis 
methods. 
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1.2 Concepts of Texture Synthesis 

Despite its importance, a precise definition of 
texture does not exist. Texture is often considered to be 
composed of a set of primitives and their spatial 
organization. More important, texture usually possesses 
an invariance property. We will use this invariance 
property as one definition of texture in this thesis. An 
observer should detect no visual difference between one 
windowed portion of a textured region and another. Thus 
texture is also a function of window size. If a 
difference over a region is detected then either the 
texture is not homogeneous (see section 2.1) or a larger 
or smaller window should be used. Windowing is very 
important when gathering statistics to be used for texture 
discrimination or texture synthesis. 

The approach to texture synthesis used in this thesis 
is outlined in Fig. 1.1. As a first step in the synthesis 
process, statistics are calculated from measurements taken 
on a parent sample texture. The statistics are then used 
to estimate model parameters. In the final step, these 
model parameter estimates are used to generate a texture 
synthesis. 

All of the digital images in this thesis are 512 by 
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512 pixels. They have either 256 gray levels (continuous 
tone) or 2 gray levels (binary) . The original parent 
texture images in this thesis have been cnosen from an 
album by Brodatz [lj. High quality prints obtained from 
the photographer were scanned and digitized at the 'JSC 
Image Processing Institute. 

The performance measure for any texture synthesis 
method is purely visual. Thus the rating of any synthesis 
system must be made by an observer and is subject to human 
variability. These facts must be remembered when 
considering the results of synthesis techniques discussed 
in this thesis and in other works. For this reason, the 
visual analysis of results in this thesis is left 
primarily to the reader. Nevertheless, general guidelines 
and trade-offs involved in texture synthesis have been 
developed in the course of this work. 

The success of any texture synthesis as well as any 
image processing technique (enhancement, restoration, 
etc.) also depends on the display medium used for final 
results. It is unfortunate that the product of so much 
work is subjected to the imperfections of recording and 
printing processes. We have attempted to minimize these 
degradations as much as possible. 
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1.3 The Stochastic Model 


The approach used to synthesize textures in this 
thesis is probabilistic (statistical) rather than 
structural because the structural approach would probably 
require a different processing method for each texture. 
Textures which are irregular and highly random are often 
difficult to analyze using a structural approach as they 
do not have "structure” as such. On the other hand, such 
textures are easily analyzed using a probabilistic 
approach. Textures which are not composed of well-defined 
non-varying primitives usually do not lend themselves to a 
structural approach. Highly regular textures which are 
usually best explained using a structural approach may 
also be studied using a statistical approach. In this 
case, the statistical model must be constructed to explain 
regular and periodic events with less importance placed on 
random elements of the model. Naturally, the structural 
and probabilistic approaches overlap considerably. 

In our statistical approach to texture analysis and 
synthesis we will use various stochastic models. Textures 
are analyzed as series of pixels to which a model is fit. 
The time or space series of successive pixels which forms 
our sample parent texture is regarded as a sample 
realization from an infinite population of such textures. 
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A model is derived which attempts to explain the numerical 
sequence of pixels in the observed parent texture. This 
stochastic model may then be used to simulate a similar 
sequence which becomes our texture synthesis. 

The power of the stochastic models presented in this 
thesis is that it is easy to use the model in a mode which 
synthesizes textures given the necessary model parameters. 
In this sense, the stochastic approach is sufficient to 
capture everything about a texture. The quality of the 
synthesis is also a measure of the amount of information 
contained in the model. 

1.4 Organization and Contributions of the Thesis 

In the next chapter, a one-dimensional texture 
synthesis model is mathematically developed. It is 
analyzed as a Markov model where the state of the system 
is considered to be defined by the sequence of previous 
pixel values. Similar texture synthesis has also been 
done by Julesz [8], Purks and Richards [3], Conners [7], 
Abend [6], and Gagalowicz [4,5], With the simple model 
presented in Chapter 2, textures with controlled 
statistical properties are generated and used to study 
human perception of texture. 

In Chapter 3, the one-dimensional model is extended 
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to tw. imensions and is used to synthesize natural binary 
textures. This work was presented earlier by Garber [9] 
and similar work was done by Lu and Fu [10]. This 
two-dimensional probabilistic model which generates 
textures based on higher-order probability densities is 
then reduced carefully to a linear autoregressive model. 
Results of both methods applied to a wide variety of 
textures are shown. 

In Chapter 4, a method to reconstruct higher-order 
densities for texture synthesis from second-order 
measurements is presented. The results confirm the value 
of second-order statistics in texture analysis. 

In Chapter 5, the autoregressive model is applied to 
continuous-tone texture synthesis. Much simpler synthesis 
models were used by McCormick and Jayaramamurthy [2] and 
Tou, Kao and Cheng [11]. The success of their method was 
not established due to the limited number of textures 
involved in their studies, however their results suggested 
that time series models could be used to generate some 
natural textures. The large, two-dimensional, linear 
model presented in Chapter 5 is then extended to a 
quadratic autoregressive form and synthesis results with 
both stationary and non-stationary noise are presented. 
The generated textures show that these models define 
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valuable synthesis methods. 

In Chapter 5, texture synthesis methods wnich 
incorporate more than one autoregressive model are 
presented. The multiple-model methods methods could be 
valuable in synthesizing textures which have very coarse 
structure and textures which are composed of subtextures. 

In Chapter 7, a texture synthesis method which 
simulates complete probability density information for use 


in the synthesis 

process is 

developed. 

This method 

is 

computationally 

burdensome 

but also 

yields the 

best 

synthesis results 

. A method 

similar 

to this will 

be 


valuable as processor speed increases in the future. 

In Chapter 3, methods for adding and removing 
non-homogeneities in texture mean and variance are 
presented . 

In Chapter 9, the measures used to estimate 
parameters in the autoregressive model are applied to the 
problem of texture discrimination. The methods differ 
from those presented earlier by Deguchi and 
Morishita [12], Kaiser [13], Pratt, Faugeras and 
Gagalowicz [14], and Haralick ejt a 1. [15]. The chapter 
illustrates two approaches to use statistics from a 
synthesis model to discriminate textures. 
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The experimental results of this thesis are largely 
visual therefore particular attention should be devoted to 
the figures. A casual reader could read section 2.1, 
Chapter 3, Chapter 5, and Chapter 7 and still understand 
much of the basic concepts of the work as well as major 
results. Chapter 10 contains a final summary and 
comparison of the models and their corresponding results. 

As a whole, this thesis presents results in texture 
simulation using methods developed herein or only briefly 
mentioned in other previous studies. The texture 
syntheses are exceptional in some cases and certainly 
notable in others. This work should encourage additional 
research in the field of texture synthesis. 

1.5 Notation 

Any variables which have different meanings within 
this thesis are clearly defined in the places where they 
are used. The term "complex" means complicated rather 
than a variable with real and imaginary parts. The terms 
"normal distribution" and "normally distributed" are to be 
taken in a statistical probability density function 
(ie.Gaussian distribution) sense. 
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CHAPTER 2 


ONE-DIMENSIONAL BINARY TEXTURE MODEL 

2.1 Introduction 

Texture is a complex image attribute that has been 
the subject of much research and is difficult to define 
precisely. The relationship between discrimination of 
textures by human observers and the mathematical 
attributes of textures has also been extensively 
researched. Models for computer discrimination have been 
proposed based on statistical parameters considered in 
some aspects to be primary texture measures. 

The terminologies in a portion of previous texture 
work have often been vague at best. As a result, the 
terms second-order and third-order have been seriously 
twisted and misinterpreted from study to study. In this 
and following chapters, we will attempt to suppress this 
confusion by carefully defining the various terms. 

The stochastic approach toward texture analysis 
considers texture fields as samples of two-dimensional 
stochastic fields. Assuming that we are dealing with 
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sampled and quantized imagery, let I(n^,n^ 2 ) denote the 
random field. Here n^ and n.^ are integers representing 
coordinates of points in the image plane. Let be the 
vector having coordinates n^ and n ^ 2 ( i .e.^= (n^^ ,n^ 2 ) ) • 
Second-order statistics are given by the set of 
second-order joint density functions 



( W 


( 2 . 1 ) 


for all possible vectors n^ and nj, where and Vj are 
the values of the random variables I(n^) and I(rij), 
respectively. In most texture work and in all of the work 
in this thesis (except for the work in Chapter 7 and 
Chapter 8 ) the random field is assumed to be homogeneous, 
that is, all orders of probability densities are invariant 
through translations. Thus, 


P->- -*■ = P-+ (? ?) 

n. ,n . n.+c,n.+c 

l ] I D 

where c is an arbitrary vector constant. As an example, 


p(v 1 ,v 2 ) = P(v 3 ,v 4 ) 


(2.3) 


where V^, V 2 , V 3 and V 4 are as shown in Figure 2.1. In 
most of our work, dummy values of random variables 
(denoted for example by V^) will be used to label pixels 
at vector location n; . 
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Given the assumptions that a texture field is 

homogeneous, the joint density functions P for all vector 

r 

separations r = n^-nj represent the most complete set of 
second-order statistics possible. The statistical 
expectation of any functions of these joint density 
functions are called second-order statistics. If a pixel 
is connected to any of its neighbors on the same row, that 
is, if we consider neighbors immediately to the left or 
right (such as Vg and Vg in Figure 2.1), then their joint 
density is called a second-order nearest-neighbor joint 
density and any statistical expectations of the joint 
density are second-order nearest-neighbor statistics. 
Nearest-neighbor densities and nearest-neighbor statistics 
are very important in this chapter as the textures to be 
generated are primarily one-dimensional. 

Similarly, third-order statistics are given by the 
set of third-order density functions 


P-+ -*■ -*• (V.,V.,V.) 

n i' n j' n k 1 1 k 


(2.4) 


Assuming homogeneity of tha texture, then 


P-V -+ -V = P-> /OCX 

n. ,n. ,n n.+c,n.-t-c,n +c U.b) 

X J K X J Jv 

for all rti and an arbitrary vector constant <?. As an 
example, 
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in Figure 2.2. The statistical expectations of any 
function of these third-order densities are called 
third-order statistics. All second-order statistics may 
be derived from third-order joint densities. In this 
chapter, third-order statistics involving adjacent pixels 
along an image row, such as the pixels, V^, Vg, Vg, in 
Figure 2.2, will be called third-order nearest-neighbor 
statistics. 

Julesz [8] created computer generated patterns with 
controlled high-order statistical properties. A 
conclusion, often referred to as the Julesz conjecture, 
drawn from his work is that texture fields differing only 
in third- and higher-order statistics cannot be 
discriminated by a human viewer. Pollack [16] showed 
later that textures whose first- and second-order 
nearest-neighbor probabilities are equal may be 
discriminated by varying the third-order nearest-neighbor 
probabilities. Purks and Richards [3] extended this 
concept to create texture patterns that differ only in 
their statistics for four adjacent points. This study 
indicates that such textures can also be easily 
discriminated. However, as was pointed out by Pratt [14], 
the second-order probability densities of the two fields 









are not constrained to be equal for arbitrary pixel pairs 
along an image line. Thus there is still some question as 
to the relationships between measured mathematical 
parameters and human d iscriminability. Later work by 
Gagalowicz [5] seems to indicate that carefully generated 
binary patterns whose second-order probability pairs are 
equal for arbitrary distances can be visually 
discriminated by human observers and therefore presents a 
valid contradiction to the Julesz conjecture. Controlling 
different statistics of a texture is often a painfully 
difficult process and as a result, most of these textures 
are generated using approaches unlike the Markov approach 
of this chapter. Often blocks of pixels are generated 
with certain properties or patterns with special 
orientation and separation are used to study the effect of 
statistical changes on the human visual system. The 
mathematical relationships between the joint density 
functions of any texture are indeed complex. 

For these reasons, we begin with one rather simple 
method of generating Markov one-dimensional binary 
textures. In later chapters, these ideas will be modified 
and extended to generate and simulate two-dimensional 
textures. 

We have studied in detail the mathematical 
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pa rameters 


involved 




relationships of parameters involved in binary 
computer-generated one-dimensional texture patterns. 
Texture patterns in this paper have been generated using 
the mathematical relationships derived herein. Methods 
have been developed to control texture statistics for both 
nearest-neighbor and non-nearest-neighbor cases. Examples 
of both types of textures are presented. Using these 
methods, numerous counter examples to Julesz's conjecture 
may be generated and are illustrated in this Chapter. 


Throughout this chapter, the P (V^,V 2 ,...,v N ) will 
denote the nearest-neighbor. Nth-order joint densities of 
our one-dimensional texture and the pixels v i» v 2'*** ,v N 
will be adjacent to one another in the sequence as shown 
in Fig. 2.3 unless otherwise stated. 

2.2 Generation Procedure 

One-dimensional binary textures represent the 
simplest form of texture possible. It is believed that 
such binary patterns force human observers to utilize 
primitive visual mechanisms in discrimination. They are 
not designed to replace or imitate natural textures but 
are experimentally valuable in deriving concepts 


concerning texture attributes due to their mathematical 
simplicity. 








Figure 2.1 Second-order 
Statistics 


Figure 2.2 Third-order 
Statistics 



Figure 2.3 One-dimensional N-grams 
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In this experiment, binary one-dimensional sequences 
with carefully-controlled transition probabilities, 
dependent on the previous four points, were generated. 
Each sequence was then broken up into shorter 512-pixel 
strings which were stacked to form the two-dimensional 
512x512 array which served as the texture pattern as is 
illustrated in Figure 2.4. Thus, the derived statistics 
are only controlled in one dimension but tne final texture 
is two-dimensional. 


We define the a priori probability of a binary 
sequence of length N by P(V^,V 2 ,...,V^) where each V^, 
i = 1,...,N is either 0 or 1. In our experiment this 
binary sequence is determined by generation parameters. 
This set of parameters, G Q (V V•••' V N ^• each of which 
represents the probability of generating a 0 after the 
contiguous binary sequence V^,V^,.•.,V N , defines the 
Markov process used to generate the texture pattern. It 
follows that the probability of generating a 1 after the 
sequence V ,V ,-' V N is 1-G o^ V l ,V 2' -' V N ) • That is ' 


vw* 


'V 


= 1-G X (V f V 2 , 


V ) 
' 1ST 


Illustration of this commonly-used texture generation 
method is given by Purks and Richards [3]. However, it 
should be pointed out that their generation parameters 
were in many cases constrained to provide equal N-gram 
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Figure 2.4 Forming a Two-Dimensional Image from 















statistics, P (V^ , V 2 , ... ) . The term "N-gram" refers to 
the Nth-order joint density function of a texture and was 
used by Purks and Richards in their study of texture. The 
term "N-gram"is frequently used in this thesis as a 
substitution for the longer terminology "Nth-order joint 
density function." A texture procedure, more general than 
that proposed by Purks and Richards, is detailed here. 

These generation parameters are actually a special 
set of conditional probabilities. Only this specific 
group of conditional probabilities is used to generate the 
texture. 


2.3 Analytical Analysis of the Markov Texture Process 


By examining the mathematics of the Markov process we 
hope to generate patterns according to a set of given 
probabilities P(V^, ,...,) which may be named the 
N-gram statistics of a specific pattern. We must 
therefore deal with the relationships that exist between 
these N-gram statistics and their generation parameters 
denoted by G(V^,•.•, V N ). Examining these relationships 
and also those between N-gram statistics of different 
lengths (that is the relationships between 


P ( V l' V 2'‘* *» v n ) and P( V i, V 2 ,•••' v n 2 ^ for a11 N 1 and N 2* 
leads us to an understanding of the probabilistic system 


involved and thereby a method of generating desired 
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texture patterns. 


In generating random texture sequences, it is useful 

to first look at a simple analogy to the process from 

which basic concepts and conclusions can be drawn. We 

might regard this generation to be equivalent to the 

experiment consisting of tossing a "smart" coin that has a 

finite memory. In this case, G q^ V i' V 2'’'*' V n^ niight 

represent the probability of tossing a "heads" given the 

previous sequence of N tosses was V,,V„,...,V . The 

12 N 

resulting string of "l’"s and "0"'s (0 is the random 

variable denoting heads, 1 denotes tails) recorded from 

this experiment is our "texture." We realize immediately 
that the texture is "determined" by this set of generation 
parameters G Q (VV^,...,V N ). Using the concept of 
conditional probability where P(A/B) is the probability of 
A given B we notice that 

G 0 (V 1 ,V 2" ' * ,V N } = P(0/V 1 ,V 2 ,...,V N ) . (2.8) 

Perhaps the most important concept derived from these 
generation parameters is that of the finite memory of the 

system. As is indicated by the notation Gq(V^, V 2. V N>' 

the probability of generating a zero depends on the string 
of binary values ' v 2 ' • • • ' v n and not those "preceding" 
• It is thereby suggested that our system has an N-gram 
memory and we will define such a system as 
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N-gram-dimensional. For example, returning to our coin 
tossing experiment, if we are in a four-gram-dimensional 
system, the probability of tossing a head depends on the 
four previous tosses only and all these conditional 
probabilities are determined by the sixteen parameters 

G 0 (V 1' V 2' V 3'V‘ 

With these concepts in mind we can find these 
generation parameters G Q (V V 2 , .. • , V N ) given the desired 
probabilities P(V^,V 2 ,...,V N ) and vice versa. The 
approach taken by Purks and Richards [3] in finding these 
N-gram statistics is based on sampling the generated 
textures. This may be seen by examining the entries in 

t 

Table I of Ref. [3]. The entries correspond to the number 
of each N-gram counted in the texture generated and the 
accuracy of such probabilities depends on the law of large 
numbers. So the true probabilities P (V ^ , . . . , V^) are 
only approximated by the output textures and this 
approximation is poor when the physical size of the 
textures is small. These estimates have greater variance 
when the true probabi1ities are small. Therefore it is 
desirable to compute the exact probabilities given the 
generation parameters of the system. 

Before proceeding further it is useful to prove the 
identity 
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(2.9) 


P(V 1' V 2.Vl } = 


P(v r v 2 ,... »V N _ 1 ,0)+P(V 1 ,V 2 ,... ,V N _ 1 ,1) 


Proof: 


p(v 1 ,v 2 ,...»v N _ 1 ,i) = 


P(V 1# V 2 ,.••#V N _ 1 )*G X (V 1 ,V 2 ,...,V N _ 1 ) (2.10) 


p(v 1 ,v 2 ,...#v N „ 1# o) = 


P — **^0 ^1 '^2 ' ’ * * '^N-l^ 

P(V 1 ,V 2 , ...,V N _ 1 )*(1-G 1 (V 1 ,V 2 ,. ..,V N _ 1 )) 


therefore 


P(V 1' V 2'''* ,V N-1 ,0)+P(V 1' V 2'*'* ,V N-1' 1) = 

P <V V 2.Vi ) * 

As a result we have the following three sets of 
equalities 

p(v 1 ,v 2 ,...,v N _ 1 ,0) = 

P(V r V 2 , . . . ») *Gg ^1 '^2 ' ’ * * ’ ^N-l ^ (2.11) 

p(v 1 ,v 2 ,...,v N _ 1 ,i) = 

P (V 1 ,V 2 , . . . fVl ).<l-G 0 (V 1 ,V 2 .V N _ X ) ) 


p(v r v 2 ,...,v N _ 1 ) = 


( 2 . 12 ) 


P (V^ ,V 2 t • • • 


G 0 (V 1' V 2'*••^N-l 5 = p (V 1 ,V 2 ,...,V N _ 1 ,0)/ (2.13) 

(P(v lf v 2 ,...,v N _ 1 ,0)+P(v 1 ,v 2 ,...,v N _ 1 ,i)). 

Equation (2.13) results from Eqs. (2.11) and (2.12) and is 
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essentially a statement of Bayes theorem for our problem. 

2.4 Texture Statistics From Generation Parameters 

Let us then consider the problem of obtaining the 
generation parameters (G's) from the N-grams (P's). 
Immediately we come to the conclusion that this is a 
trivial problem. One might merely use Eq. (2.13) to 
deduce the generation parameters. However, the equations 
derived so far do not indicate the extensive relationships 
which exist between the N-grams. 

Equation (2.13) is not invertible so it is not useful 

in obtaining the N-gram statistics, P(V 1 ,...,V N ) of a 

sequence given the generation parameters. As was stated 

above, once the generation parameters are defined, a 

texture may be generated using those parameters and the 

N-gram statistics are determined. We also know that once 

a complete set of N-gram statistics P(V..,...,V ) are 

N 1 

defined for some N , the N-gram statistics P (V ,...,V ) 

1 y 1 N 2 

may be resolved using Eq. (2.12) for all n 2 <N i* Given the 
generation parameters of a system, G^(V^,V 2 ,...,), we 
can analytically determine the N-gram statistics, 
P (Vi, V 2 ' * * * » V n ) resu lting texture. 

The solution to this problem of finding N-gram 
statistics given generation parameters may be found by 
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considering the generation procedure as a discrete state 

Markov process. This approach is readily seen wnen 

considering the generation parameters G rt (V. , V_ , . . . , V ) as 

0 12 N 

transition probabilities. If we consider a 

two-dimensional system with P (0,0),P(0,1),P(1,0) and 

P(l,l) and generation parameters G(0,0), 0(0,1), G(1,0) 

and G(l,l) we may define our system as composed of four 

possible states (0,0), (0,1), (1,0) and (1,1). If the 

system is in state i at the Kth observation and in state j 

at the (K+l)th observation then we say that the system has 

made a transition from state i to state j at the Kth stage 

of the generation process. In our example an observation 

is taken at each generation of a single new binary value 

and the state is determined by the values of the last two 

binary numbers generated. As an example, consider the 

sequence 0,1,1,0,0. We might consider the system to be in 

the (0,1) state at the start which may represent the Kth 

stage of our generation process then a transition is made 

to the (1,1) state at the (K+l) stage. These transition 

probabilities are determined by the generation parameters 

of the system. We also note that our N-dimensional system 
N 

has 2 possible states. As the transitions from each of 
these 2 N possible states to each of the 2 N possible states 
is fixed by our generation parameters we may form a 
transition matrix T whose elements t(i,j) represent the 






probability of a transition from the ith state to the jth 
state. If T is the transition matrix of a regular Markov 
chain, then there is a unique probability vector p which 
has positive coordinates and satisfies 

T T p = p (2.14) 

This same vector p may be computed by taking any row of 
the matrix 

T q (2.15) 

as q approaches infinity [17]. The vector p represents 
the vector of steady state probabi1ities. In our case it 
contains the desired probabilities P(V^,V 2 ,... ,) . 

It is important to realize that this theorem holds 
for regular Markov processes. If there is an integer q 
such that every element of the matrix in Eq. (2.15) is 
stricly positive then the process is regular. Some 
processes are not regular such as absorbing Markov chains 
[17]. When any element of the transition matrix is equal 
to one along the diagonal, the process is said to be 
absorbing given that the system may begin in any state. 
This could happen if Gg(0,0) = 1 for example (a series of 
0's would be generated in this case). For the purposes of 
our discussion we will assume that 
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for all V-^. This is a sufficient but not necessary 
condition for the process to be regular. 


Applying 

these 

concepts to 

a two- 

d imensional 

system 

we obtain the 

transition matrix 







Final 

State 




- 

00 

01 

10 

11 



00 

o 

o 

o 

o 

1-G q (0,0) 

0 

0 


Initial 





State 

01 

0 

0 

G 0 (0 ,1 ) 

1-G q (0,1) 

(2.17) 


10 

O 

r—H 

o 

o 

1-G 0 (1,0) 

0 

0 



.11 

0 

0 

G 0 (l,l) 

1-Gq (1 ,1 )_ 



The first row contains the transition probabilities 
from state (0,0) to states (0,0), (0,1), (1,0) and (1,1) 
in that order. The following set of equations results 
when Eq. (2.14) and Eq. (2.17) are combined: 


G Q (0,0)-1 

0 

Gq (1,0) 

0 


P(0,0) 


0 

1-G q (0,0) 

-1 

1-G 0 (1,0) 

0 

X 

P(0,1) 


0 

0 

G Q (0,1) 

-1 

G 0 (l,l) 


P(1,0) 


0 

0 

1-G q (0,1) 

0 

-God,!). 


P(U) 


_ 0 _ 


As the above system is singular , we may form an 








equivalent non-singular set by replacing any equation with 


P(0,0)+P(0,1)+P(1,0)+P(1,1) = 1 (2.19) 

by using the fact that p is a probability vector. Solving 
this system gives the desired two-gram statistics 

p(v L ,v 2 ). 

Examining these generation parameters further we find 
that the same N-gram statistics may be generated by 
generation parameters of a different dimension. For 
example, consider the generation parameters in Table 2.1. 


TABLE 2.1. EQUIVALENT GENERATION PARAMETERS 


Set 1 


Set 

_3 

G q (0,0) = 0.05 

O 

o 

o 

o 

= 0.05 

Gq (1,0,0) = 0.05 

G q (0,1) = 0.07 

g 0 (0,0,1) 

= 0.07 

Gq(1,0,1) = 0.07 

G q (1,0) = 0.92 

G 0 (0,1,0) 

= 0.92 

G (1,1,0) = 0.92 

Gq (1,1) = 0.75 

Gq(0,1,1) 

= 0.75 

Gq(1,1,1) = 0.75 


Notice that the probability of generating a zero 

following a v i» v 2' v 3 ^ oes not depend on V^. The values, | 

| 

G 0 ^ V 1 ' V 2 ' V 3 ^ secon d set indicate that the system is- } 

memoryless beyond two previous generation steps. We may 

write 
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G 0 (V 1 ,V 2' V 3 ) = P(0/ Vl ,V 2 ,V 3 ) 


( 2 . 20 ) 


p ( o/v 2 ,v 3 ) = g q (v 2 ,v 3 ) 
for all V 2 and 

It follows that, according to Eq. (2.11), the 
P(V^,V ,...,V ) are also determined in our example for N>2 
given P(V^,V 2 ) and G Q ( V ]_» V 2 )* Thus we have 


p(v 1 ,v 2 ,v 3 ) = p(v 1 ,v 2 )*g v (v 1# v 2 ) 


( 2 . 21 ) 


P(V 1' V 2' V 3 ,V 4 ) = P(V r V 2 ,V 3 )*G v (V 1> V 2 ,V 3 ) 

= P(V 1 ,V 2 ,V )*G V (V 2 ,V ) 

4 


( 2 . 22 ) 


We conclude that given any set of generation parameters 

G 0 (Vi,V 2 .V N ) we may form a set of generation 

parameters Gg (V-]_, V 2 , . . . , V N ) , M greater than or equal to N, 
according to the rule 


G 0 (V 1 ,V 2' ' • ’ ,V M-N' V M-N+1'* * • 'V " 

g o ( Vn + i.V 


(2.23) 


that generate an equivalent set of N-gram statistics and 
therefore equivalent textures. 


When desiring to generate textures according to a 
given set of N-gram statistics, P (V^,V 2 ,...,V ) we must 
realize the set of constraints imposed on the set. For 
example, 
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(2.24) 


r 




•V 


1 


for all N. Returning to the set of equations used to 
determine the 2-gram statistics in matrix form we realize 
that, by adding the first two rows and last two rows of 
Eq. (2.18), P(0,1) = P (1,0). In fact, by considering the 
set of equations arising from the set of equations derives 
from the generation systems of higher dimensions we find 


P(V X 


,v 2 ,v 2 ,v 2 ,. 
p(v 3 ,v 2 ,v 2 


•"W = 


,v 2 ,...,v 2 ,v 1 ) 


(2.25) 


This implies that many constraints 
N-gram statistics. For example, in 
system containing P(V^,V 2 ,V 3 ) 

P(0,0,1) = P (1,0,0) 
P(0,1,1) = P(1,1,0) 


are present on the 
the 3-gram-dimensional 


(2.26) 


but also by Eq. (2.12) and the fact that P(0,1) = P(1,0), 


P(0,1,0)+P(0,1,1) = P(1,0,0)+P(1,0,1) 


(2.27) 


and 


P(0,0,1)+P(1,0,1) = P(0,1,0)+P(1,1,0) 


By definition we also know that 
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P(0,v 1 ,v 2 ,-..,v N _ 1 )*G 0 (0 / v 1 ,v 2 ,...,v N _ 1 ) 


+P(1,V 1 ,V 2 ,...,V n _ 1 )*G 0 (1/V 1 ,V 2 ,.../V N _ X ) (2.28) 
= P(V 1 ,V 2 ,... ,v N _ x ,0) 


and 


p(0,v 1 ,v 2 , 

+P(1 / V 1 ,V 
= P(V lf V 


2 

2 


' V N-1 ) * [ 1 " G 0 ( 0 ' V 1' V 2'- 

••' V N-1 ) * [ 1 _G 0 ( 1 ' V 1' V 2 


'Vl )] 


• ’ ,V N-1 )1 


(2.29) 


Combining Eqs. (2.28) and (2.29) with Eq. (2.12) 

P(0 ' V 1' V 2'*‘• ,V N-1 )+P(1,V 1' V 2'‘'* ,V N-1 ) = 

p(v lf ...,v N _ x ) 

Equation (2.30) holds for all N>2. 


(2.30) 


Four constraining equations exist when Eqs. (2.24), 

(2.25) and (2.30) are reduced to orthogonal form for this 

3-gram-dimensional system. This implies that we have four 

degrees of freedom when choosing the eight values 

P ^1 ,V 2 ,V 3^ along with the obvious constraint that 

P(V ^,^2 > v 3 )>0• This is precisely the number of degrees of 

freedom in the set of g q^ V 1 ' V 2 ^ which have only the range 

constraint of Eq. (2.16). For higher N-grain-dimensional 

N-l 

systems, P (V-^, .. . 'V always has 2 constraints on it 

from Eqs. (2.12), (2.24), (2.25) and (2.30). Thus we see 
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that Eq. (2.13) holds in a degrees of freedom sense. 

In conclusion, a method of determining N-gram 
statistics from generation parameters using the concept of 
a Markov chain was developed and a set of linear equations 
describing the constraints on the N-grams was presented. 
Using this approach to generating texture patterns a large 
variety of textures may be easily generated and examined 
using a minimum amount of effort. 

2.5 Texture Moments 

The above equalities and inequalities provide a full 

understanding of the texture generation system in 

probabilistic terms. Still further conclusions can be 

derived from them. From the above we see more clearly 

that the generation parameters G(Vj,V 2 ,...,V N ) determine 

the texture completely and thus define the N-gram 

statistics P(V^,...,V M ) for all M. Also for a given set 

of P(V.,...,V ) there can exist an infinite number of 
i M 

generation parameters, G (V^,V 2 ,...,V ), which would 
generate many textures with such statistics if N>M-1. 
Provided the constraints on the statistics P(V^,...,V ) 
are met just one texture could be generated if N = M-l or 
perhaps none at all if N<M-1. Thus textures with equal 
first, second, third and fourth nearest-neighbor 
probabilities can be generated. 
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Parameters thought to be useful in texture 
discrimination may also be easily developed. For example, 
joint moments about the mean defined as 


TTY T 

E[(x 1 -p 1 ) 1 (x 2 -u 2 ) 2 (x 3 -y 3 ) 3 ...(x k -, k ) k ] (2.31) 

y 

where ^ r^ is the order of moment [18]. The rth moment 
of x^ is defined as 

E ( x i> xJf( Xl .x k ) (2.32) 

X 1 X 2 X k 

where f is the joint probability distribution of the x^. 
From our binary textures we could define the following 
parameters: 

U - E{x 0 > P ' x i> - 

0-P(0)+l-P(l) = P (1) 
a 2 = E{ (Xq~u ) 2 } =^~^(x-y) 2 f (x) = 

(O-P(l)) 2 P(0) + (1-P(1)) 2 P(1) 

= P(l)-P(l) 2 

E { (x Q -y) 3 } =^(x-y) 3 f(x) = (2.33) 

(O-P(l) ) 3 P(0) + (1-P(1) ) 3 P(1) 

= 2P {1) 3 -3P(1) 2 +P(1) 

E{(x Q -y)(Xj-y)} _ P(l,l)-P(l) 2 

a — 2 — ■ 5 

a P(l)-P(ir 
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1 


P(1,1,1)-3P(1,1) *P (1)+4P(1) 3 
(P(l)-P(l) 2 ) 3/2 

where P(l), P(l,l), P(l,l,l) represent nearest-neighbor 
(N-gram) statistics although this can be changed to 
include non-nearest-neighbor statistics thus creating new 
texture parameters. The above parameters are useful in 
discrimination therefore only when textures differ in 
their (3-gram) or shorter statistics. 

2.6 Constraining Second-Order Statistics 

We describe now a method which allows 
non-nearest-neighbor statistics to be controlled using the 
relationships developed in the Sections 2.3 and 2.4. 
Because second-order probabilities are of interest we 
investigate the conditions required to assure equality of 
second-order statistics for non-nearest-neighbor 
statistics. If we denote G g (V V 2 » • • • > v j^) as our 
generation parameters and P (V V 2 , • • • , V N ) to be their 
associated N-gram statistics then fo'r the second-order 
statistics of one texture to be equal to another for the 
(N-l)st neighbor distance, 


E{ (Xq-Vj) (x 1 ~u) (x 2 ~u) } 

3 

a 


3 







(2.34) 


I] •••L P a <V l' V 2 . V M ,V M+1. V 


v 2 V N -1 (2-34) 

E - ’ ' E Pb(V l' V 2"-- ,V M+1' • • " V N } 

' V 2 V N-1 

for V^,V N e{0,l} where P a , represents the N-gram 

statistics for the first and second texture respectively. 
It can also be shown that 


E ’* E P(V 1' V 2. V M. V = 


V 2 V N-1 


E-E p (v i ,v 2" • •' 


v 


(2.35) 


V 2 V N-1 


V, 


where 


kXi‘ V< ' 1) k ■ Wm-Vm+i .\-l” 




Vj Vj =0 

Recall that the N-gram statistics, P, are a function 
of the generating parameters G. If the second-order 
statistics are to be equal for two textures regardless of 
neighbor distance then Eq. (2.35) must hold for all N. It 
was previously believed that combining Eq. (2.34) and 
Eq. (2.35) yielded a non-redundant non-linear set of 


equatio 

ns that 

would imply 

that 
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having 

equal 
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statistics 

must 

have 
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parameters [19 
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distributions cannot exist. He then proceeded to work 
with textures of four grey levels. 


In the next section, we will show that binary 
textures with equal second-order statistics for all 
distances can be generated using a carefully-chosen set of 
generation parameters and that these texture contadict the 
Julesz conjecture. 

2.7 Experimental Results 

We may define second-order statistics for 
non-nearest-neighbors as 


(2.36) 


p< w -Z •••X) Z •••E p<v i' v 2. v j. 


V 


V. 


V j-1 Vl V N 


Purks and Richards [3] attempted to demonstrate that 
textures equal in second-order distribution could be 
generated with visually detectable differences. However, 
they merely held second- or third-order statistics equal 
between the two textures for small j, while varying 
second- and third-order statistics for longer j, that is, 
they allowed the second-ordor statistics for 
non-nearest-neighbors to be unequal over some distance. 
Examples of this type of manipulation are shown in 
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Figs. 2.5-2.10. The corresponding N-grams are shown in 
Table 2.2 and Table 2.3. The first column of each 

subtable contains the analytic N-grams for texture (a), 
top, and texture (b), bottom, which are found by solving 
Eq. (2.14) and Eq. (2.24) given the generation parameters 
in the third column. The middle column contains the 
actual N-gram statistics as measured from the generated 
textures themselves. The input, analytically-solved 

parameters will not be equal to the output statistics as 
the generation process is random. In other words, the 
output statistics are based on measurements taken on a 
sample from a population of textures having the 

characteristics defined by the input parameters. 

Figures 2.5 and 2.6 show pairs of statistics having 

equal first- and second-order nearest-neighbor statistics. 

That is, P (V ,V ) = P. (V,,V„). There is also an internal 
a 1 2 b 1 2 

equality for these textures which may be expressed as 

p a (Vf, v 2) = Pb (V l' V 2 ) = 1/4 for a11 V 1 and v 2 that are 

nearest-neighbors. Visual differences are apparent 

between the pairs. Figures 2.7 and 2.8 have texture pairs 
which have equal third-order, nearest-neighbor statistics 
both within and between the pairs, that is, 

P a (V 1 ,V 2 ,V 3 ) = P b ( v i»V 2 ,V 3 ) = 0.125. The bottom texture 
in Fig. 2.7 is a coin flip, purely random texture with the 
probability of both black and white equal. It is 
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Figure 2.7 One-dimensional Figure 2.8 One-dimensional 
Texture Texture 
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Table 2.2 N-GRAMS AMD GENERATION PARAMETERS 
FOR FIGURES 2.5 - 2.8 
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Table 2.3 N-GRAMS AMD GENERATION PARAMETERS 
FOR FIGURES 2.9 - 2.12 
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interesting to note that this texture is visually quite 
similar to the bottom texture of Fig. 2.6, yet the 
generation parameters and N-grams differ significantly. 
Thus two textures with differing N-gram statistics may be 
generated which are visually similar. Figure 2.9 shows a 
pair of textures which . have equal fourth-order, 
nearest-neighbor statistics both within and between the 
pairs. That is, 


P a (V l ,V 2' V 3' V 4 ) = P b (V l' V 2' V 3' V 4 ) = °' 0625 • (2 ' 37) 


The two are easily discriminated. Figure 2.10 shows two 
textures with between-equa 1 , fourth-order nearest-neighbor 
statistics. The N-grams are not equal within however. 


Figures 2.11 and 2.12 have texture pairs which are 
similar in many ways to the pair in Fig. 2.6. Both are 
counter examples to Julesz's conjecture that the eye is 
sensitive to only first-order and second-order probability 
distributions. Each has a texture pair where second-order 
statistics are equal for all distances between the texture 
pair. Precisely stated, p a ( v i» v j) = p b^ v l ,v j^ f° r a33 3* 
These textures are generated using a set of restrictions 
discussed in Appendix A. 


The textures in these figures are pseudo-randomly- 
generated. Each texture was tested using a 
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goodness-of-fit procedure to insure that textures with the 
desired N-grams were generated. This was done for 1-grams 
to 10-grams. The procedure is outlined in Appendix B. 
The goodness-of-fit depends on the pseudo-random number 
generator used, as a poor generator can be guaranteed to 
yield poor results in this type of experiment. The 
pseudo-random number generator used is detailed in 
Appendix C. 

2.8 Conclusions 

The one-dimensional binary patterns generated give 
rise to some basic concepts concerning textures and their 
discrimination. First of all they indicate the use of 
moments and similar statistics is not optimal at least in 
the nearest-neighbor sense as many textures have equal 
moments but are visually quite different. However, it 
should be pointed our that this may only be characteristic 
of some artificial textures and that moments could serve 
as good discrimination parameters in many real-world 
applications. Secondly, the results indicate a close 
relationship between second-order non-nearest-neighbor 
statistics and human discrimination. The counter-examples 
to the Julesz conjecture indicate that N-grams and 
higher-order statistics may be valuable in identifying and 
discriminating some textures. Use of these texture 
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measures depends on factors such as discrimination 
accuracy desired, cost factors for statistics measure and 
the nature of the textures involved. 


In later chapters, investigation of two-dimensional 
textures will be pursued as these correspond more closely 
to natural scene textures. The texture generation task 
will be approached from a simulation rather than a pure 
synthesis point of view. The complexity of controlling 
two-dimensional statistics in a synthesis process and 
attempts to study their effects on human discrimination 
and interpretation of textures is a problem which requires 
careful analysis and is beyond the scope of this work. It 
is very possible that no applicable, clear-cut results 
could be obtained from such a study. However, by 
simulating textures, models and processes for generating 
similar-looking textures are derived which may be useful 
in other applications and some simple statements 
concerning human discrimination of textures can also be 
made. 
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CHAPTER 3 


TWO-DIMENSIONAL BINARY TEXTURE MODELS 

3.1 Introduction 

In this chapter, the concepts used to generate 
one-dimensional binary textures are extended to the 
two-dimensional case. This model is then used to simulate 
natural binary textures. A method for choosing the pixels 
in a non-contiguous generation kernel based on a linear 
model is described. This is an important concept in much 
of the work presented in this thesis. The method for 
collecting N-gram statistics is discussed and practical 
problems arising in this process are investigated. 
Results of natural texture simulation using this method 
are presented. Finally, the linear model which was used 
to determine the kernel pixels of greatest value in the 
N-gram generation process is used to generate binary 
simulations. This will lead us to the application of the 
linear model to continuous-tone textures in Chapter 5. 

3.2 The 2-D Binary Markov Model 

In the investigation of natural phenomenon once a 
researcher collects enough data he tries to imagine a 
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process :ich accounts for the results. The construction 
and development of a mathematical model is often the best 
way to do this. done. In some cases, the model may be 
extraordinarily complex, in others, exceedingly simple. 
In most cases model acceptance cannot be based on "truth" 
as the true generating phenomenon is too complex or simply 
unknown and so it is based upon model usefulness and "how 
well it works". Such "working" models for texture are 
presented in this chapter because they have the ability to 
simulate some natural textures. 

If one can synthesize and simulate natural textures 
adequately by using some proposed model the criterion of 
usefulness and workability for that model is met. A 
researcher may then also apply the model to problems of 
texture identification and discrimination with 
justification. With any set of texture measures required 
for simulation, the information content of the measures is 
viaually indicated by the quality of the simulation and 
not merely by the percent of correct versus incorrect 
classifications when those textures are applied to the 
discrimination problem. By adding features important to 
texture simulation better methods of texture 
discrimination and identification can possibly be found. 

Many early texture studies involved the use of binary 
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textures generated by one-dimensional Markov processes. 


Such work was 
one-dimensional 
generated line 
pa rameters 


presented in Chapter 2. 
models a large vector of 
by line using a set of 


In these 
pixels was 
generation 


where 


V (V 1' V 2.V 

N+l 


G V ( V 1/V 2 ,...,V N ) = P(V N+1 /V 1 ,V 2 ,...,V N ) (3.1) 
N+l 

and P(A/B) represents the probability of A given B. In 
the above notation each represents a generated pixel 
which has value 0 (black) or 1 (white). Each pixel value, 
then, depends on the N pixels previous to it. A 
two-dimensional texture image is then formed by breaking 
up the large vector of pixels into shorter strings and 
stacking them one on top of the other (see Fig. 2.4). 
This procedure for large images nearly insured image row 
independence (unless N was large) thereby creating only 
horizontally oriented textures totally unsuitable for 
simulating natural two-dimensional textures. 


By allowing N to increase exceeding the short string 
line length, two-dimensional (vertical and horizontal) 
dependence may be induced into the generating process. A 
pixel value then depends not only on the pixels previous 










to it on the same line but also on the pixels above it 

(see Fig. 3.1(b)). Thus, textures could be generated as a 

time sequence in television raster scan fashion. In 

theory, texture dependence could be extended _ad infinitum , 

however practical considerations concerning the actual 
• N 

generation process show us that 2 generation parameters 
must be accounted for. As a possible solution to the 
storage problem we can choose to ignore, all but N of the 
previous pixels in our generation process and we can allow 
the pattern of the 's to become flexible. This idea 
will be discussed in later sections of this chapter. 

Throughout the remainder of this thesis, the set of 
pixels, VVs, on which the next pixel, V N+1 , depends will 
be referred to as the "kernel" of the synthesis process. 
The pixel V N+1 will be referred to as the "eye" of the 
kernel. 

In order to estimate p < v i' V 2 ''* ,V n } for a fixed 
pattern V ^,V 2 ,...,V , M substrings (samples) of length 

N are taken from a parent substring of length M+N-l and 
the number of occurrences of the specific pattern 
V ,V ,...,V are counted, then divided by M. This is 
equivalent to estimating the probability density function 
of a random variable by the histogram of a set of samples. 
Ignoring boundary conditions, linear unbiased estimates 
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(P's) of the P's may be defined as 


p ( v i/ v 2 


'V 


1 

M 



N 

n 6 (I (k+j ) -v, ) 

k=l K 


(3.2) 


where 


6 (V v V k ) 


0 

1 


if V j * V k 


if 


V . 
D 


V 


k 


and I(i) represents the ith element of the one-dimensional 

texture string from which the parameters are to be 

estimated. Equation (3.2) assumes that the V. are 

i 

contiguously located in order along a line. 


This idea of estimating N-grams, P(V.,V ....,V ) from 

12 N 

a sample parent texture may be extended to the 
two-dimensional case. A histogram of occurrences of each 
pattern of (V^, , ...,V^) is made by passing the 
two-dimensional kernel in Fig. 3.1(b) over the 
two-dimensional sample parent image. The tally is then 
divided by the total number of sample patterns observed to 
obtain P(V^,...,V ). As was stated earlier, 

two-dimensional synthesis is merely an extension of the 
one-dimensional case ignoring boundary conditions of the 
two-dimensional image. 


Although it was not explicitly stated in Chapter 2, 
the generation parameters of a texture may be estimated 
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These 


for any given set of N V^'s from a parent texture, 
statistics have the property 


E [G 


V (v i» v 2 ' 

N+l 


,V )] = G (V ,V ,...,V M ) 
N V N+1 1 2 N 


where 


V < V 1- V 2.V ■ P( V V 2. V N-Vl )/ 

N+l 

(p(v 1 ,v 2 ,...,v n ,o)+p(v 1 ,v 2 ,...,v N>1 )) . 


(3.3) 


3.3 Seeding the Generation Process 

In the one-dimensional texture pattern generation of 
Chapter 2, a long binary vector of pixels is broken up 
into short strings which are stacked on top of each other 
(see Fig. 2.4). When synthesizing this one-dimensional 
string, a seed of four binary values is actually required 
to begin the process for a 4-gram-dimensional system. 
When the Markov process is regular and non-absorbing with 
none of the generation parameters equal to one or zero the 
process rapidly reaches a steady state independent of the 
seed. 


A similar thing happens in the two-dimensional case 
but here the seed is larger and forms a two-dimensional 
frame around the synthesized image texture as shown in 
Fig. 3.2. With a two-dimensional texture generation 
kernel such as that shown in Fig. 3.1(b), the synthesized 
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(a) One-dimensional (b) Two-dimensional 

Figure 3.1 Texture Synthesis Kernel 



Figure 3.2 Synthesized Image Texture and 
Seed Region 
























image texture may be generated without using the bottom of 
the frame so the next pixel at each generation step 
depends only on pixels above it. 

The seeding process may be handled in a variety of 
ways. The simplest approach would be to randomly generate 
the seed once for the whole image. In this case the pixel 
values in the seed frame of Fig. 3.2 remain the same 
throughout the generation process. A second approach 
would require the random generation of each pixel, V., in 
the generation kernel that fell outside the synthesized 
image texture region at each step during the generation 
process. Using this method, the pixel values in the seed 
frame change at each pixel generation. A third method 
would involve wrapping the image around such that the left 
edge of the synthesized image texture joins with the right 
side as in Fig. 3.3. In this case a random seed is 
required only to begin the process of the top of an image. 
A final method involves the use of another texture, 
usually a previously generated texture or the parent 
texture, as part of the seed, rather than noise. This 
method reduces the noise which otherwise occurs around the 
edges of the synthesized image texture. 

Regardless of the seeding process, all texture 
synthesis methods developed in this thesis normally 
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converge to a steady state within 5 to 20 pixels of the 
border of the image. This was confirmed by repeated 
studies of convergence effects on texture simulations. In 
most cases, this narrow region is not noticeable and is 
included as a part of the result. In some critical 
applications these edges could be thrown away. 

3.4 Kernel Selection Using The Linear Model 

We will refer to the V.'s on which the next pixel, 
V N+ ]_, depends as the kernel of the generation process. 
Geometrically speaking, the V^'s form a kernel "shape" or 
"pattern" which may or may not be spatially contiguous. 
For example, in Fig. 3.4 a generating kernel shape is 
shown where the V 5 pixel directly depends on only some of 
the pixels in its surrounding neighborhood. In this case, 
V,_ may be generated based on the values of pixels 
Vf, v , V 3 and V 4 but is directly dependent on no other 
pixels in the neighborhood. This does not imply that 
is not related or correlated with its other neighbors. In 
fact, the relationships between V^, V 2 » V^, V^, and 
will determine other interrelationships. 

A non-contiguous neighborhood of V^'s is used as it 
allows a more parsimonious model for texture generation to 
be chosen. An analogy is in simple linear regression (as 
defined by Draper[20]) where independent variables which 
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do not contribute to the prediction or estimation of the 


dependent variable are dropped. In texture generation 

this allows the model to be estimated by fewer parameters 

and makes the generation-synthesis process more efficient 

by reducing the number of computations required. When 

generating textures based on N-grams, reducing N reduces 

N 

the amount of storage required for 2 generation 

parameters. By allowing the kernel of 's on which V N+ ^ 

depends to be non-contiguous, the range of dependence in a 
distance sense is increased over that which would be 
allowed with a contiguous kernel containing the same 
number of V^'s. This is very important to obtain the 
larger structure apparent in many textures. Reducing the 
number of pixels in the model also relieves us from the 
complex numerical problems of inverting matrices of 
unwieldy size, a necessary step in linear model parameter 
estimation discussed later in this chapter. We would, for 
example, not expect our V N+ ^ pixel to depend on a pixel 
where the spatial separation between V N+ ^ and is large. 
If that distance is small, however, we would expect a 
large dependence. 

The method for choosing the proper independent 
variables (V^'s) to be included in the generation process 
requires special attention. We wish to choose the best 
subset of N variables from a larger finite neighborhood of 






T variables, where N<T. Evaluating such subsets and their 

corresponding models requires a criterion. Texture 

results for each possible model could be visually examined 

and compared and the V^'s of the model corresponding to 

the visually most pleasing result could be chosen. 

T 

However, ( N ) model evaluations must be done using this 
approach. For a simple search through T = 40 points with 
N = 12, 5.5 billion models would have to be evaluated! 

This approach is therefore impractical and so a 
sub-optimal approach which yields a good but not 
necessarily the best set of 's for our model must be 
used. 


If we view this problem as one of predicting a 
dependent variable, Vi f rom a large set of independent 
variables, V^'s, then the standard linear model approaches 
may be applied. In a statistical sense, independent 
variables are values that can be observed but not 
controlled and dependent variables are affected by changes 
in the independent variables. Thus the value of dependent 
variables is said to depend on values or changes in 
independent variables. The linear model is just one 
approach to explaining the relationship between 
independent and dependent variables. 

In most linear model applications, the criterion used 
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to evaluate the fit of the model to data is mean-square 

error. It is desirable in most cases to choose a model 

which minimizes the mean-square error. The problem in our 

case is to choose a subset of \A ' s of size N from a set of 

1 s of size T such that the linear model employing those 

's produces the minimum mean-square error when compared 

to all other possible models containing N V^'s. This 

T 

cannot be done without examining all ( N ) models again 
however, suboptimal approaches producing a very good, but 
not necessarily the best fit are available. One method 
employing a forward selection procedure is described in 
the following sections. 

As a final note it should be observed that we are 
choosing the kernel for a non-linear N-gram-based on the 


fit of a linear 

model 

to sample data. 

This approach 

is 

admittedly ad 

hoc 

and is chosen 

for simplicity 

and 

computational ease. 

Regardless, it is 

believed that 

the 


kernel chosen by this approach is very good if not the 
best. If the value of a particular is important to the 
prediction of it will have a high partial correlation 
coefficient when it is examined for entry into the model 
during the forward selection proceedure. This is true in 
all cases when the pixels are binary-valued and is usually 
true when the pixels are continuously valued. Still, as 
we are not considering all (£) possible sets of pixels in 








the neighborhood, the best model will not always be found 
[ 20 ] . 


The linear model which we use to determine the V^'s 
in the kernel may be expressed in linear regression form 
as 

Y k = xj E k k = 1,2, . . . ,M 

where 


Y k = V N + l,k 


1 

V l,k 

V 2 ,k 


*-V, 


(3.4) 


N,k 


B is an (N+l)xl vector of unobservable parameters and 
is an unobservable random variable such that E[e k ] = 0. 
The sample number is denoted by k and there are a total of 
M samples. We can also define matrices X and Y as 


X = 


- - 


— 

+T 

x i 


Y 1 

X 2 


Y 2 

• 

Y = 


• 


• 

T 


• 

. V 


_ y m. 


(3.5) 


-*• T -1 T 

The most common estimator of 8 is (X X) X Y, the 


least-squares estimator. so given a parent texture wnich 
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we desire to simulate, an estimation of the texture model 
parameter 3,3, may be obtained yielding the generation 
model 

Y = xt (3.6) 

It is very important to note that the estimation of the 
amount of noise present in a texture is usually obtained 
by measuring the amount of variation in the parent texture 
which is unexplained by our model. This is expressed 
numerically in the mean-square error. Therefore, as in 
any modelling process, any shortcomings or inaccuracies of 
the model will appear to be "noise" (unexplained variance) 
and hence the amount of noise will increase. 

3.5 Correlation and Partial Correlation 

We may define the mean and variance of a variable as 

M v = E [V^] (3.7) 

and 

% = Et(v r‘ , v 1 )2 ' <3 - 8) 

Similarly we may define the covariance of two variables 

and V as 
2 

Y„ „ = E[(V,-y v )(V,-v v )] 

V 1 V 2 1 1 ^ 2 

(3.9) 


* 'I'Vz' - Vv ., 1 
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And their correlation coefficient as 


'V, V 


12 


V 1 V 2 


a V a V 
V 1 2 


(3.10) 


Using the above definitions we may then define the partial 
correlation coefficient of variables V-^ and V 2 after both 
have been adjusted for as 


w 


V. 


PV 1 V 2~ (PV 1 V 3 )(PV 2 V 3’ 

^ V 1 V 3 2 ^' P V2 V 3 


(3.11) 


Each of the above parameters has a corresponding statistic 
(estimate of a parameter of a population given an 
observable sample of that population) which is chosen to 
meet some desirable set of a criteria such as sufficiency, 
consistency and unbias of estimate under certain 
conditions. Given that the ^'s are samples from an ith 
order multivariate normal distribution the maximum 
likelihood estimates for the above parameters are 



(3.12) 




<v i,k-V 2 


M 


(3.13) 
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I>i,lT V l ) (V 2,iTV 


ViV2 " {[E (v i.^i )2 ][Z (v 2,rt>] 2 } 1/2 


(3.14) 


r — (r ) (r ) 

V 1 V 2 V 1 V 3 V 2 V 3 


' V lV V 3 


Vi-4 v 


(3.15) 


" V 1 V 3 


V 2 V 3 


Second-order partial autocorrelation may be found in a 
similar manner using 


SvVXvVVV^ 


V 1 V 2- V 3 V 4 




V lV V 3 




(3.16) 


v.v..V, 

2 4 3 


Higher-order partial correlations may be found by 
extension of the above [21]. 


3.6 The Forward Selection Procedure 


One approach to finding a good subset of independent 
variables is known as the forward selection procedure. 
This method inserts N variables into the model 
one-at-a-time. The order of insertion is determined by 
using the partial-correlation coefficient as a measure of 
the importance of variables not yet in the equation. The 
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basic p idure is as follows. First we select the V i 
most correlated with V N+1 denoted as , and we find the 

1 A 

first-order linear regression equation v N+ i = a i v i +a 2• 
We next find the partial correlation coefficient of \T 
(for all j / i) and V N+ ^ (after allowance for ). 


Mathematically this is equivalent to finding the 
correlation between the residuals from the regression 

A 

V N+ ^ = a^V^ +a 2 and the residuals from another regression 

V. = b V. +tv (which we have not actually performed). 
] i ^ 

The Vj with the highest partial correlation coefficient 

with v N+ i» v ^ * now selected and a second equation 

V = c V. +c V. +c_ is fitted to the data. This 
N+l 1 2 ±2 ^ 

process continues. After V. ,V. ,...,V. are in the 

1 2 l q 

regression the partial correlation coefficients are the 
correlations between the residuals from the regression 

A 

V N+1 = f(V^ ,V. ,...,V^ ) and the residuals from a 

1 2 q 

regression V i = f. (V i ,V ^ , . .. ,V i ) (j^q). The V i with 

j 3 1 2 X q X j 

the highest partial correlation coefficient is now 

selected for entry into the linear model. The process is 

continued until N 's are entered into the model. 


The final N variables chosen by a forward selection 


procedure are not guaranteed to be the optimal set but 
given the logistics of the selection procedure, the 
solution obtained is usually close to optimal. 








r 


One o£ the most common procedures for implementing 

the forward selection process numerically utilizes 

Doolittle decomposition [22], The Doolittle decomposition 

may be used to find the inverse of the correlation matrix, 

R , and the estimates of linear model coefficients as each 
P 

variable is entered in the model. The correlation matrix 
merely consists of the set of correlation coefficients 
r v y as defined by Eq. (3.14). It is then factored into 

i j 

the product of two triangular matrices 


R 

P 


L U 
P P 


(3.17) 


where L is lower triangular and U is upper triangular 
with ones on the diagonal. Partial autocorrelation 
coefficients may be obtained easily from elements of this 
matrix during its decomposition at each step. For further 
details and examples on the forward selection procedure of 
the decomposition process see Beyer [18] and Draper and 
Smith [20]. 


3.7 Statistics Collection 


In practice, N, the number of pixels in the 
generation kernel, is often chosen by computational limits 
imposed by finite processing capability or finite computer 
memory. The idea of parsimony would also urge the 
selection of the smallest N possible. In the most general 
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texture synthesis (stochastic) model which utilizes 

N 

N-grams we find that as many as g storage locations are 
required in order to collect the data needed to synthesize 
a texture from a parent where g is the number of grey 
levels in the texture. This approach calls for N to be 
less than 17 for g = 2 (a binary image) if we have only 
2 ^ = 131,072 storage locations in memory. Approaches to 
"stretch" this limitation have been investigated and will 
be discussed in a latter section. If we have 8 grey 
levels then N must be less than or equal to 5 as 

5 

3 = 32,768. Thus, in synthesizing textures using an 
N-gram approach, processor storage capability is the major 
limiting factor. 

Determining which pixels will be included in the 
generation kernel requires an estimate of the linear model 
defined in Eq. (3.4). To do this the X and Y matrices of 
Eq. (3.5) may be used or the correlation matrix of the 
kernel points must be estimated using a parent texture. 


The 

elements 

o f 

the 

*k 

vector of 

Eq. 

(3.4) are 

obtained 

by passing 

the 

kernel window over 

the sample 

parent texture and 

recording 

the 

pixel value 

corresponding 

to the 

pos i t ion 

o f 

each V ^ 

i n 

the window. 

The 

kernel is 


assumed to be completely within the boundary of the parent 
texture. A sample may be taken at all possible positions 






in the parent texture or a random sample of points may be 
chosen if the parent texture is very large to reduce the 
number of computations required. In actual practice, no X 
matrix is ever formed. To obtain the correlation matrix 
of the sample set, R, (and from that 6 ) we need only the 
surn of squares cross products, and V ^, over the sample 
set according to Eq. (3.14). 

The elements of the correlation matrix for many 
kernel patterns often contain redundancies. For example 
the spatial relationship between the pair V-, and V 2 and 
the pair V 3 and V 5 in Fig. 3.4 is the same. Estimating 
the correlation for these two pairs of points from a 
sample will yield nearly equal values if the sample size 
is large and the overlap of samples used to estimate the 
correlation values is large. In short, 

r V.V 9 = r V-.V. (3.18) 
12 3 4 

In more mathematical terms, let Ifn^n^ denote the 
random texture field where n ^ and n^ are integers 
representing the coordinates of points in our sampled and 
quantized image. Let n be the vector of coordinates 
( n i, n 2 ) • From Section 2.1, second-order or 2-gram 
statistics are given by the set of joint density functions 

P -> (3.19) 
n ,m J 
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for all possible vectors n and in, where and Vj are the 
pixel values of the random variables I (n) and I (in) , 
respectively. If the random field is homogeneous, that 
is, invariant through translations then 


P = P 

-y -v 

n,m n-m, 0 


from Eq. (2.2). 


(3.20) 


If we assume that our two-dimensional texture is 
homogeneous and stationary we might propose that 
Eq. (3.20) is true by definition. A word of caution is 

necessary at this point. Estimating the elements of the 
covariance or correlation matrix using the assumption of 
stationarity and homogeneity can yield correlation 
matrices having negative determinants, in violation of the 
fact that non-singular covariance matrices must be 
symmetric positive definite. The violation occurs because 
the sample is not homogeneous and stationary. Therefore 
the type of assumption expressed in Eq. (3.20) should only 
be used for estimation when sample sizes are very large 
and are known to be homogeneous and stationary. 

The assumption of Eq. (3.20) can be very powerful in 
a computational sense for large samples as the calculation 

A 

of covariance (or correlation) matrices can be time 
consuming. The complete covariance matrix Cor a 
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contiguous 10 x 10 (100 element) kernel contains 10,000 


entries. Of these, 4050 are redundant by simple symmetry 


as / j M V = / j V V which implies 
l j 3 l 


r = r 
V • V • V • V • 
i j j i 


Computing all cross products and sums to estimate this 
matrix from a 512 x 512 image requires over 1.025 billion 
multiplications and 1.05 billion additions. Utilizing the 
concept of homogeneity, this may be reduced to 0.05 
billion multipiications and 0.075 billion additions which 
is significantly lower. 


Once the covariance matrix of a kernel is determined, 
the linear model containing the V^'s of the kernel may be 
obtained. It is advisable to make the kernel large as 
more texture information is contained over large 
distances. In many cases, there are relationships between 
pixels separated by great distances especially if the 
texture is coarse or highly regular (periodic). In 
practice, however, large matrices may be numerically 
ill-conditioned during a decomposition or inversion 
process and so they must be avoided. Also computer 
storage limitations must be considered. In this study, no 
matrix larger than 100x100 was decomposed or inverted for 
these reasons. This constraint required a multipi e-pass 
approach. 


First, 100 V^'s (usually closest to the kernel eye. 
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V ,) were chosen for examination. They were entered into 

the linear model in a forward selection manner until it 

was determined that entry of additional V^'s would be 

insignificant. Insignificance is indicated both 

statistically by considering the reduction in the sum of 

squares due to the entry of a variable into the model [20] 

and also by the value of 6 which approaches zero so the 

i 

variable V. becomes insiynificant. On the second pass, 
variables not tried in the first pass are examined and 
tested for possible entry into the model. Again, those 
which are significant are kept and the others are 
discarded. On the next pass, any variables not examined 
in the first two passes may be examined. The process 
continues until all V^'s have been tested for possible 
entry into the linear model at least once. 

This multiple-pass process, as the forward selection 
procedure which it employs, is not guaranteed to produce 
the best model but will provide an excellent model 
nevertheless. The final model will also be parsimonious. 
This will reduce the number of parameters in our model and 
aid in the efficiency of the generation process. 

3.8 Results 

Once the points for the kernel are chosen based on 
the linear model derived using the methods described in 
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the previous sections, estimates of the generation 

parameters tor the texture are obtained using concepts 

discussed in section 3.2. Practical considerations 

require us to limit N, the number of pixels in the kernel, 

to 12 to 13 depending on the processor storage available 
N 

as 2 values must be stored. These G's are then used to 
generate each pixel along a row, row by row until a 
complete two-dimensional texture is obtained. For eacn 
pixel the appropriate generation parameter estimate is 
found and a uniformly-distributed pseudo-random variable 
is generated. Based on these two values, a black pixel 
(0) or white pixel (1) is generated. 


In practice, not all of the generation parameters may 
be estimated when N is large because all possible patterns 

of v i» V 2 ' * * * ' v n' V N+1 ma ^ not * 3e P resent if the sample 
image or there may be few of them. Smaller samples can 
cause inaccurate estimation of the G's as the variance of 
our estimate is larger and therefore the expected error of 
our estimate is larger than would be expected with a 
larger sample size. In these cases it is important to sum 
over the least significant kernel elements and estimate 

by G„ (V.,V.,,,...,V„). In our 


°V <V V 2.V 

N+l 

study, this was done if the sample size to compute 


-v N+1 ( Wi.V- 


G (V ,V ,...,V ) was less than 10. The variable i is 
V N+1 1 2 N 

increased until this condition is met. 
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'i-.. uure simulations using this method are shown in 
Figs. 3.5(b)-3.14(b). Visually, the results are very 
good. As the estimated texture generation parameters are 
approximated using statistics gathered from the full 
parent texture, non-homogeneity in the parent texture will 
cause an "average" texture to be synthesized. The 
simulation of straw (Fig. 3.7) is poor because of its 
non-homogeneous nature (specifically the directionality of 
the stalks in different parts of the image) and exhibition 
of detail (specifically individual non-conforming single 
stalks). A similar observation may be made with respect 
to the parent textures of grass and water but in these 
textures the non-homogeneity is not so pronounced. As we 
are attempting to synthesize textures and not merely 
"image code" the parent textures, details and 
non-homogeneities will be lost in the synthesis process. 

The bark texture is among the most difficult to 
simulate due to its very unusual macro-structure. Still, 
the N-gram simulation looks remarkably similar to the 
original when windowed regions 20 to 40 pixels square are 
obse rved. 

* 

The kernels used to generate these N-gram simulations 
are shown in Table 3.1. A clear-cut relationship between 
the textures and the kernels found using the linear model 
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(a) Original Texture 


(b) N-gram Simu 



(c) Linear Model Simulation 













































(a) Original Texture 


(b) N-gram Simulation 



(c) Linear Model Simulation 
Figure 3.9 Leather 
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Original Texture 


(b) N-gram Simulation 



(c) Linear Model Simulation 


Figure 3.12 Wood 
















Figure 3.13 Pigskin 
























Table 3.1 TWO-DIMENSIONAL N-GRAM TEXTURE 
SIMULATION KERNELS 


GRASS BARK 
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approach is not always evident. Some kernels, such as 

those for wool and water, seem to have a vertical or 

horizontal structure possibly resulting from the structure 
of their corresponding textures while others, such as the 
kernels for wood and raffia, do not. Still, the 

aon-contiguous kernels produce excellent results. 

3.9 Extension of .M-Grarn Model 

The N-gram model may be extended beyond the fixed w 

set by limited processor storage. The requirement that G 

must be estimated based on no fewer than q (which was set 

to 10 in our study) samples already reduces the storage to 

a maximum of M/q non-red undent parameters where M is the 

total number of samples in the parent texture image. For 

5 

a 512x512 image, M is approximately 2.5x10 . Given 

q = 10, this implies that a maximum of 25,000 generation 
parameters estimates must be stored. 3y increasing q, 
further reduction is possible. Meanwhile N is permitted 
to increase without bound to sample size limitations. For 
frequently occurring patterns, larger N's allow increasing 
dependency over larger distances. This is most desirable 
in regular and coarse-structured textures. 

.Searching for the proper generation parameter at each 
step in this type of process is complex. In the earlier 
N-gram storage, each pattern of 0's and l's actually forms 


x 
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a binary address to the location in memory of the desired 
G. In the extended case a sort, search, or a series of 
comparisons along with some intelligent preprocessing is 
required. Efficiency is reduced. 

To illustrate the effect of model extension the 
texture raffia was used. Figure 3.14(a) shows the 
original and Fig. 3.14(b) is the synthesis obtained using 
a kernel containing N = 14 pixels. Figure 3.14(c) is the 
linear model synthesis. Figures 3.15 (a) ,(b),(c) were 
obtained using three different texture kernels with N = 22 


points. 

Far 

more structure 

in 

these extended 

model 

versions 

is 

apparent. This 

is 

expected as at 

each 


generation step, the next pixel is allowed to depend on 
pixels further from it. As the pixels in the kernel 
become more widely spaced the synthesis becomes more 
structured but small, local regions often become more 
distorted and less raffia-looking because the information 
used in the synthesis process is more global than local. 

3.10 Linear Model Generation of Binary Textures 

The process of choosing the V ' s to be present in the 

l 

texture generation kernel described in earlier sections of 
this chapter actually yields a simple linear model which 
can also be used to generate binary textures. The model 
which results from the determination of the generation 
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kernel may be expressed in equation form as 


V N + l,k " 8 l V l.k + 8 2 V 2,k + --- +e N V N,k + V £ k 


(3.21) 


or more simply as 


V N + 1 ' 8 lV 8 2 V 2 + --- + W V e ’ 


(3.22) 


Once the estimates of the 3^'s are known, a pixel V N+1 may 
be calculated from a set of given values plus an error 
e. In one-dimensional analysis this is sometimes known as 
the autoregressive time series model [24J. For binary 
a value of V N+ ^ will be produced which is non-binary. To 
generate binary data using this model will therefore 
require quantization. 


In the N-gram approach to texture simulation, the 
randomness of the texture is induced by the generation of 
a uniformly-distributed pseudo-random variable during the 
generation process. The comparison of this value with the 

A 

estimate of the generation parameter, G, yields the next 
binary pixel. A similar type of randomness must occur in 
the generation of binary textures using the linear model 
of Eq. (3.22). This randomness is expressed in the model 
in the error term e. 


We can obtain an estimate of the distribution of e in 
the same manner as we estimate the g's of the model. This 
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may be done by applying the model to the sample data from 
which it was derived and observing the errors. That is, 
the linear model kernel is passed over the parent texture 

A 

image and at each point a V is calculated based on 

N+l 

A 

Eq. (3.22) without the error term. Then V -V is 

N+l N+l 

calculated where V is the actual value of V in the 

N+l N+l 

parent texture. The histogram of the values can be used 
to estimate the distribution of £. As one step further we 
could assume that has some known distribution such as 
Gaussian or normal, and merely estimate the parameters 
necessary to define this distribution. In the normal 
distribution case, only the standard deviation (or 

variance) of needs to be estimated. The mean of e is 

zero in the linear model, least-squares distribution. 

Our generation process then consists of the 

calculation of ^ ^ V.; + B Q to which we add a random, 

normally-distributed error term e and this value is then 

quantized to 0 or 1 based on comparison with 0.5. Results 

using this generation method are shown in Figs. 3.5(c) 

through 3.14(c). In these figures, N was allowed to be as 

N 

large as 70 as only N coefficients (not 2 ) need to be 

stored along with , the estimate of the error standard 

deviation. 

The kernels used in the linear model simulations are 
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p 


shown in Table 3.2. The linear model simulations are 
slightly inferior to the N-gram simulations but the 
degradation is far less than we would expect from such a 
massive compression of information (which is approximately 
2 to 70). The results were good enough to encourage the 
application of the linear model to continuous-tone 
textures. 
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Table 3.2 TWO-DIMENSIONAL BINARY LINEAR MODEL 
TEXTURE GENERATION KERNELS 


GRASS 




BARK 

cP a 
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CHAPTER 4 


ALGEBRAIC RECONSTRUCTION TECHNIQUE MODEL 

4.1 Introduction 

Julesz's conjecture [8] that second-order statistics 
are sufficient for human visual texture discrimination 
provides a useful estimate of the amount of information 
necessary to reconstruct a texture field. Although 
counterexamples to that conjecture have recently been 
found [4,25], and are shown in Chapter 2, it is a good 
first-order approximation. Examples of the use of that 
upper bound for texture analysis can be found in [14,26]. 
Therefore it is very tempting to use it for synthesizing 
natural texture fields. That is, we may attempt to 
simulate textures based on generation parameters estimated 
from 2-gram statistics over various distances. This 
approach requires fewer statistics to be collected from a 
parent texture as only 2-grams versus N-grams are 
collected. 

This chapter illustrates that we must "invent" 
higher-order statistics to use the Markov generation 
coefficient approach for texture synthesis if we limit our 
knowledge of the original random field to second-order 
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statis ts. We will demonstrate how this can be done 
using Algebraic Reconstruction Techniques (ART). 

4.2 Problem Definition 

The stochastic approach toward texture analysis 
considers texture fields as samples of 2-D stochastic 
fields. Second-order statistics are given by the set of 
second-order joint density functions (2-grams), P(V^,Vj). 
As in Chapter 3, we also assume for simplicity that the 
random field is homogeneous, that is, invariant through 
translations. Given the set of joint density functions, 
(P(V if Vj), for all spatial relationships between and 
Vj, our problem is to synthesize a texture field with the 
same second-order statistics. 

There are many ways of carrying out such a synthesis. 


We use 

a 

television 

raster 

type of scanning (left 

to 

right, 

top 

to bottom) 

and 

the generation kernel 

of 

Fig. 3. 

1(b) . 

Even if we 

assume 

finite memory, namely that 


intensity at point V N+ j depends only on intensities at 
points located in some finite neighborhood, we see that 
second-order statistics are not sufficient to generate the 
process. Indeed, assuming a memory of order N, the 
intensity at V is computed using the conditional 
probability or generation coefficient 

Gy (V^jVji • • • iVjj) = ‘ ’' (4.1) 
N+1 
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which involves (N+l)th-order joint density functions. 
Therefore if N is larger than 1 and we are given only 


second-order statistics, P(\T,\T), we have to "invent" 
higher-order densities P (V^,V 2 ,...,). Mathematically, 
the problem can then be stated as follows: given N random 
variables V ]_'*‘ - ' V n such that for every pair 
l^i,j<N and i f j, we know the joint density function 
P(Vi,Vj), find a function P(V^,...,V N ) which satisfies 


P ( V 1 , V 2 » • • • / V N ) > 0 for all V ]_'**** V N 


E-E E-E E-E p<v i.v- 


v, v. . v.,, v. . v. , v. T 

1 l-l l+l ]-l 3+1 N 


* v j »••• » V N )- 

(4.3) 


P(V i ,V j ) . 

Here the P(V ,V )* s are sometimes called the marginals of 
i 3 

P(V ,...,V ). Assuming quantization with g levels, the g N 
IN 

unknowns P (V , . . . ,V ) can be stacked as a vector and 

I N 

conditions (4.2) and (4.3) correspond to a linear 
programming problem: 


I Ap = rn (4.4) 

p > 0 (4.5) 

where vector is obtained from the functions P(V^,V^) and 
matrix A of size (( 2 )g^)xg N contains only ones and zeroes. 
For any reasonable values of N and g, this is a set of 

linear equations and inequalities of fairly large 
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dimension and the usual solution techniques such as the 
simplex method [4] become very limited. 

4.3 Solution Through Algebraic Reconstruction 

Techniques (ART) 

The ART algorithm was introduced by Gordon, Bender 
and Herman [27] for solving the problem of 
three-dimensional reconstruction from .projections in 
electron microscopy and radiology. This is a 
deconvolution problem in which a function in a 
higher-dimensional space is estimated from its 
experimentally measured projections in a lower-dimensional 
space. For an excellent review of those techniques see 
Gordon [28]. 

The problem stated in Eqs. (4.2) and (4.3) or (4.4) 
and (4.5) is precisely of this form, where the projections 
are the second-order joint density functions. ART is 
therefore directly applicable. The simple intuitive 
interpretation is that each projected density is thrown 
back across the higher-dimensional region from whence it 
came, with repeated corrections to bring each projection 
of the estimate into agreement with the corresponding 
measured projection. 

Formally, we use an iterative scheme defined by 
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P (q+1) (V, , . . . ,v„) = P (g) (V, , . . . ,v„)+tc (q) (V, , . . . ,V X1 ) 


(4.6) 


For all values of V^,...,V N , for q = 0,1, 


where the correction term c^) (V]_ , . . . , Vjj) is given by 


c ^ (V 1 ,...,V N ) - 
N-l N 




(4.7) 


( 2 J i=l j=i+l 


and P (V^, Vj ) is the actual marginal measured, for example, 

% (q) 

from an original texture field. Here P (V ,V ) is the 

i j 

marginal corresponding to the reconstructed density at 
iteration q. 


We may express this in words as follows. The 

iterative process may be started with all reconstruction 

, 'v (0) , . n 

elements set to a constant (P (V ,...,v ) =l/g for all 

1 N 

(V ,...,V ). in each iteration the sum of the differences 
1 N 

between the actual and the reconstructed marginals is 

N—2 

computed and evenly divided amongst the g 
reconstruction elements. If the correction is negative, 
it may happen that the calculated density becomes negative 
at some points. This problem can be alleviated by using a 
modified iteration scheme defined by 


p (q+1) (V r . . . ,V N ) = MAxjo,P <q) (Vj^ .. . ,V N )+tc <q) (Vj^.. .,V N )} 

Mq + 1) (4 - 8> 

therefore guaranteeing P ^ 0 (constrained ART [28)). 


It is of course necessary to determine when an iterative 
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algorithm has converged to a solution which is optimal 
according to some criterion. This is in turn related to 


the problem of finding the optimal value t of t. Various 
criteria for convergence have been devised [28], For 
simplicity, we chose the mean-square error 


£ = ||<£ 
q 


(q) 


2 

2 


i -*■ 'v (q) | , 2 


(4.9) 


between the measured and calculated marginals where ||•| L 
is the usual euclidean no/m and in is defined in Eq. (4.4). 


To derive the optimal step size, t, for each 
iteration we rewrite Eq. (4.6) in vector form as 




(4.10) 


Multiplying both sides of Eqs. (4.10) with matrix A 
(Eq. (4.4) we obtain 


£(q+D = ^ (q) 


m VM, + td 


(q) 


(4.11) 


where 


J<q) . A 3<q) 


(4.12) 


and subtracting the actual marginal vector m from both 
sides of Eq. (4.11) yields 

| |S (q + l , ||2 „ ||j(q). t a<q)||2 

= t 2 l |d <c3) [ ||-2t3 (q) -e <q) + | [e <q> | || ! 4 -“’ 


Therefore the error r at iteration o+l is minimized for 

q + 1 
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(4.i4: 


= gtq) ,a(q> 

tq " Il3 <q) l!2 

where • denotes the inner product. 

A dual approach also explored is based upon the 

analysis of Eq. (4.3) in the Fourier domain. It car 
easily be shown that the initial problem stated by 
Eqs. (4.2) and (4.3) is equivalent to an interpolation 

problem in the Fourier domain. The major drawback of this 

approach is the difficulty to ensure the positivity of the 
inverse Fourier transform of the interpolated function. 
Therefore this method was not pursued even though it may 
be the case that "good" interpolating functions will 
alleviate that problem. 

The basic philosophy of the two approaches just 

discussed is that Nth-order joint density functions are 

"invented" to satisfy exactly the constraints stated in 

Eq. (4.3). Their obvious disadvantage is the high 

dimensionality (g N for an Nth-order joint density 

function) of the data that is to be stored compared with 

N 2 

the usually lower dimensionality ((^)g for the 
second-order joint density functions) of the data that is 
effectively used. 
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4.4 Results and Conclusions 


The iterative process defined in Eq. (4.6) may be 
halted at any number of iterations, q, and a texture may 
be generated using the value of p at that point. However, 
it should be kept in mind that the success of a texture 
synthesis depends on making the error c^ as small as 
possible and that the texture generation process is 

sensitive to this error. It has also been found by 

% 

experimentation that the p contains many values which are 
set to zero by implementation of constrained ART. This 
tends to cause the Markov texture generating process to 
become absorbing, which causes patches of white and black 
or streaks and lines to be generated. This is eliminated 
by setting those values which are zero to a small non-zero 
value, 6 , in the generation process. 

Using the above concepts, texture simulations of the 
binary textures water (Fig. 4.1(a)) and raffia 
(Fig. 4.2(a)) were generated (Figs. 4.1 (b) ,4.2 (b) ) . 
Textures similar to those in Chapter 3 employing actual 
Nth-order statistics (N=14) were also generated 
(Figs. 4.1 (c) ,4.2 (c)) and are included here for reference. 
The in the generation kernels for the ART model and the 
N-gram model were chosen using the ideas described earlier 
in Chapter 3. 
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(a) Original Texture 


(b) ART-generated Simulation 



























(c) N-gram Simulation 
Figure 4.2 Raffia 





















At the onset of investigation of the Algebraic 
Reconstruction Model it was hoped that this approach would 
be useful for generating continuous-tone textures. 
However, using it to generate binary textures revealed 
that the convergence for the iterative process was very 
slow even with the optimal step size t. Each iteration 
required much computation and the storage required for the 
N-grams was large. More than eight hours of CPU time on a 
DEC KL-10 were required to execute the large number of 
iterations required for a visually pleasing solution. 
Generating textures based on estimated N-grams which is 
detailed in Chapter 3 is probably more efficient and less 
complex computationally. This work did lead to other 
texture simulation models (see Faugeras [29]). 

The results using algebraic reconstruction are nearly 
equivalent to the results using complete N-grams as 
second-order statistics collected from binary textures 
contain a great deal of information. Close examination of 
the textures generated using algebraic reconstruction on a 
high-resolution display device reveals a high random noise 
level. The computational requirements and the final noise 
level indicate that the N-gram method of generating 
textures is much less complex and yields better visual 
results than the method employing algebraic 
reconstruction. 
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CHAPTER 5 


CONTINUOUS-TONE LINEAR TEXTURE MODEL 

5.1 Introduction 

In this chapter, the concept of using a linear model 
for generating binary textures that was briefly discussed 
in Chapter 3 is extended to multi-grey-level (256-level) 
or continuous-tone textures. Pictorial results of the 
application of this model to simulation of a variety of 
textures are presented. 

The application of the linear or autoregressive model 
to time series processes has been extensive. These 
applications, which range from weather forecasting to 
stock market predictions, primarily utilize the models and 
concepts introduced by Box and Jenkins [24], Many authors 
have expanded and elaborated on their approaches 
[30-38,62]. Some researchers have applied these ideas to 
texture simulation, in fact, the Box-Jenkins 
autoregressive model is one of the very few approaches to 
texture synthesis presented thus far. 

McCormick and Jayaramamurthy [2] were perhaps the 
first to make a notable attempt to simulate natural 
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textures using this approach. Their work consisted of a 
discussion of the Box and Jenkins autoregressive (AR), 
moving average (MA) and autoregressive integrated moving 
average (ARIMA) models including estimation of model 
parameters and adequacy of model fit. (These terms are 
later defined in this chapter.) A very simple model was 
then used to simulate two very similar textures which 
closely resemble the wood texture of this study by filling 
in the holes of a parent texture using the derived model. 
Only two textures, both exhibiting a wood-grain-like 
structure, were used. Similar work was done later by Tou, 
Kao and Chang [11]. Unfortunately, the results of their 
simulation of these textures were displayed using a 
printout of Chinese cYiaracters and so the degree of 
success of their method is unclear. The appearance of 
texture synthesis results on a computer printout will 
confuse most observers unaccustomed to such crude image 
displays. The models were again very simple and contained 
no more than three terms in the linear model summation. 
Deguchi and Morishita [12] attempted to use the linear 
model to segment and partition textures. Their approach 
was only partially successful. 

In the above simulation attempts, the models used 
were simple. The process of collecting statistics and 
estimating parameters is complex. In some cases, previous 
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authors attempted to use the complex Box and Jenkins ARIMA 
model which leads to difficult model parameter estimation 
if the number of model elements is greater than two or 
three. 


In our study, the simpler autoregressive model is 
used and is allowed to contain a large number of 
parameters. This is possible using the assumption of 
homogeneity (stationarity) combined with the forward 
selection process of choosing non-contiguous generation 
kernels as described in Chapter 3. These models are 
extended further by allowing second-order autoregressive 
models and non-stationary noise. Results of texture 
simulations using these models are included in this 
chapter. 

5.2 The Linear Autoregressive Model 

In Chapter 3 the linear autoregressive model, used to 
determine the elements of the generating kernel, was 
expressed as 

Y k = *k^ + e k k = 1,—,M (5.1) 

where 

Y = V 
k N+l,k 
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and 


X k = 


1 1 k 
2,k 


L v 


N, k 


Here 3 is an (N+l)xl vector of unobservable parameters and 

e is an unobservable random variable such that E[e,] = 0. 
k k 

The sample number (index) is denoted by k and M is the 
total number of observations. We can also define the 
vectors Y and t and the matrix X by 


r ±T 


r 


- 

X 1 


Y 1 


e i 






X" 


Y_ 


E 

2 


2 


2 




-> 


• 

Y = 

* 

e = 

• 

• 


• 


• 






L X M 


- y m- 


• e M‘ 


and our model may be expressed as 

Y = X$ + e (5.3) 

In equation form, dropping the k subscript, the model 
becomes 
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Sums and sums of squares leading to the calculation of the 
correlation or covariance matrix of the parent texture are 
obtained by passing any chosen generation kernel pattern 
over the texture. From this matrix, the least-squares 
parameter estimate of 6 is obtained. The multiple-pass 
forward selection process described in Chapter 3 leads to 
a final linear autoregressive model which is then used to 
generate textures. 

5.3 Autoregressive and Moving Average Models 

In this section, we will introduce the general linear 
model as defined by Box and Jenkins [24J and Grabill [39] 
and discuss the relationships that exist between it and 
the autoregressive and moving average models. When the 
autoregressive model is extended without bound it is 
essentially equivalent to the general linear model and the 
moving average model is a subset of general linear models. 
Allowing for a large model reduces the complexity of 
parameter estimation and allows easy selection of a 
generating kernel and model. 

Many of the equations in this section also assume 
that the texture is one-dimensional and that the 
generating kernel is contiguous. That is, V N+1 follows 





V . This is consistent with notation presented earlier if 
images are expressed in lexicographic [40,41] notation as 
one-dimensional vectors. The non-contiguous kernel may be 
expressed as contiguous if the Bj_ are zero in the model. 

The output of a linear filter whose input is white 
noise may be described using the general linear time 
series model 


% 


V N+1 = £ N+1 + Vn + *1 £ N-1 + Vn-2 + '-* 


= £ N+1 + 


X^j e N-j 


j=l 


(5.5) 


where {y^ + ^=V N+ ^- y and p is the mean of the process, 

assumed to be stationary here. Thus V N+ ^ is a weighted 

sum of present and past values of the white noise process 

2 

e k * e k usually has zero mean and constant variance o £ . 
The auto-covariance is also defined as 


Y t * EU kW 



(5.6) 


Notice that Eq. (5.5) may be written in a different 

form 
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(5.7) 


^N+l * 0 V 1T 1^N-1 + TT 2^N-2 + '" + C N+1 


OD 

= ^ ; 71 • ‘ 

' 3 N-] 


+ £ 


N+l 


j=0 


Equation (5.7) may be derived from Eq. (5.5) as 


'N+l 


^N+l E-/^j £ N-j 


j=0 


N 


E^j £ N-j 


j = l 


Therefore 


^N+l = E N+1 +I ^0 


v TVe* • 

N Z —4 J N-j 

j = l 


+ Z^j E N-j 

j = l 


£ N+l +lJ; 0 V N + ZL/^j E N-j 

j = l 


Similarly 


'N-l ^N-l S^j^N-j 

j=2 


=c> ^ = e 

N+l E N+1 


+ *oV (1 '*0 )< 'l r N-l +(1 '*0 ) E*j E N-j 

j = 2 


= e N+l + *oV (li O ) *l 


^N-l E^j E N-j 

3 = 2 


+ (1- ^0 ) 2 lt ' 

3=2 


(5.8) 


(5.9) 


(5.10) 


(5.11) 


(5.12) 

j C N-j 
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oo 


j = 2 

And so, by continuing this process Eq. (5.7) is found. 

It may also be shown that Eq. (5.5) can be rewritten 

as 


V N+1 ^ {B)e N+l 

where B is the backward shift operator 


(5.13) 


Be. = e. . 
k k-1 


(5.14) 


B J e, = e 


k k-j 


(5.15) 


and 


(MB) * (5.16) 

j = 0 

where 4> = 1. ^ (B) is often called the transfer function 

of the linear filter. 

As certain constants or parameters must be estimated 
from the sample data available it is sometimes important 
to minimize the number of parameters required to 
accurately represent a process. This simplifies analysis 
of a model and reduces the number of required 
computations. The general linear model containing an 
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injiini number of terms is of little practical value. 

Therefore, if is often expedient to allow the general 

linear time series process to be reduced to a model in 

which the current state of the process may be expressed as 

a finite aggregate of previous values of the process and a 

driving error value e . This autoregressive (AR) model 

N+l 

may be written as 


a. 


K+l ~ *C)V hK-l + ^2^N-2 + ‘ * + VN-p^N-! 


(5.17) 


We may define the autoregressive operator <J>(B) as 


4><B) = (1-<|) 1 B- <I> 2 B - 


B P ) 


(5.18) 


and thus Eg. (5.17) may be rewritten as 


^ ^ V N+1 G N+1 ’ 


(5.19] 


The moving average (MA) model may be written as a 
special case of Eq. (5.5) where only the first g+1 of the 
'P weights are non-zero. This process is 

V N+1 = e N+l ” 9 0 c N _ 9 l e N-l" **' ~ 0 p e N-p 

g (5.20) 

e N+l“ Z^ 0 j G N-j 
j = 0 

As in the case of the autoregressive model, we may write 
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the moving average operator as 


6 (B) = 
and thus Eq. 


(l-e.B-e.B 2 - ... - 9 .B q -e B q+1 ) 
0 1 q-1 q 

(5.20) may be rewritten as 



0(B)e 


N+l 


(5.21) 


(5.22) 


A model incorporating the finite-term concept of both 
the autoregressive and moving average model may be written 
as 


= 4>,,V +<}>.. ,+...+ <J> ^. T +E„ 

N+l ON 1 N-l p N-p N-l 


0 O e N _ 6 l £ N-l“ 6 q E N-q 


(5.23) 


or 

<M b) V n+ i = 6(B) En+1 . (5.24) 

This mixed autoregressive-moving average (ARMA) processing 
can be thought of as the output from a linear filter 

whose transfer function is the ratio of two polynomials 
9(B) and <p (b ) when the input is white noise e . 

An Example 

To show the relationship between an autoregressive 
process and a moving process we will consider the simple 
example 
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% _ 

V N+1 ~ e N+l“ 6 0 e N 


which is a simple moving average process, 
equation indicates that 


The 


£ N e 0 E N-l 


This implies that 


C N = V e O e N-l 


Substituting Eq. (5.27) into Eq. (5.25) yields 


^N+l 0 O^N~ e O £ N-l +£ N+l 


and by a similar substitution of 


£ N-1 ^N-l +e 0 £ N-2 


into Eq. (5.28) we see that 


Vl = -'oV'SVl-<',-2Sl 


Continuing this process yields 


co 


Thus, a 
infinite 




(N + l) - j + f 'N+l 


j = l 

finite moving average process may be writte 
autoregressive process. 


(5.25) 

above 

(5.26) 

(5.27) 

(5.28) 

(5.29) 

(5.30) 

(5.31) 

as an 
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5.4 Parsimony Between Models 


The above equations show how the general linear time 
series model is related to the autoregressive, moving 
average and autoregressive-moving average models. The 
concept of parsimony suggests that is is usually desirable 
to express a model with as few terms as possible. 
However, it has been shown that if we are willing to 
sacrifice the concept of parsimony and deal with infinite 
(or large) general linear models we can use the infinite 
autoregressive model (Eq. (5.7)), also known as the 
general linear model, to describe any general linear time 
series process. 

In theory, it is possible to approximate as closely 

•» 

as desired any general linear time series process with a 
finite autoregressive process of order p by allowing p to 
increase until the desired closeness is obtained. In 
application, however, the estimation of model parameters 
tk of Eq. (5.7) may be less accurate than desired due to 
the noise in the system. That is, the noise of a 
system will cause error in the estimation of the to the 
extent that some tt^'s are believed to be zero or are 
estimated so poorly that an inaccurate model is developed. 
It is not clear to what extent the accuracy of a model is 
improved in such a case by using an autoregressive moving 







average model. It is clear, however, that going to such a 
model causes increased complexity in the parameter 
estimation process and induces difficulty into the 
hypothesis testing process which is simple in an 
autoregressive model case. Parameter estimation for the 
ARMA model requires a multiple-pass, extensive and 
computational1y-complex iterative process which is 
complicated by the required extension of our model in a 
two-dimensional image to many points and very large sample 
size. For this reason, the ARMA model was not used for 
texture simulation in this thesis. 

5.5 Results 

The linear (autoregressive) model of Eq. (5.4) was 

used to simulate 3 variety of natural textures. 

Stationary, independent Gaussian noise was used to drive 

the synthesis process. The variance of the noise was 

estimated by applying the model to the sample data and 

observing the prediction errors. Images resulting from 

the fitting of estimated models to sample data are shown 

later in this chapter. These errors, which are often 

called residuals, are pixels formed by the difference 

V -V , where V , is the actual observed pixel value of 
N+l N+l N+l 

trie sample parent image and is the corresponding 

fitted value obtained by use of the linear model. The 


112 






standard deviation of these errors can be measured and 


1 


used as the standard deviation of pseudo-random 
normally-distributed noise in the generation process. 
Actually, this information can also be obtained during the 
decomposition of the covariance matrix. 

The number of pixels in each generation kernel, N, 
varied from 30 to 60. The generation kernels used for 
each simulation are shown in Table 5.1. The simulation 
results are shown in Fig. 5.1(b) through Fig. 5.11(b). 

These simulations indicate that the linear model 
using stationary gaussian noise produces acceptable 
simulations of a variety of textures including grass, 
wool, leather, sand and water. As with the binary model, 
the simulation of bark (Fig. 5.2) shows absence of macro 
texture. The non-homogeneities present in the straw 
texture (Fig. 5.3) cause an "average" straw to be 
generated which is very similar to the binary texture 
simulation. The cloth texture (Fig. 5.4) is composed of 
two subtextures and therefore simulations made using 
statistics measured over the whole image will be a mixture 
of the two subtextures. The simulation of raffia 
(Fig. 5.11) has good structure as the linear model kernel 
is large but sharp edges present in the original texture 
are absent in the simulation. The same is true of the 
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Table 5.1 TWO-DIMENSIONAL LINEAR MODEL TEXTURE 
GENERATION KERNELS 


GRASS 


BARK 



WOOL 


m 


% 
CZZ3 
o □ 


o 

D 


cd 
, □ 


□ 


D O ^ 

CD □ 

n n O n 


rf> ° 

a a 


* 


□ 
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Figure 5.1 Grass 
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(c) Second-order Linear 
Model 


(d) Second-order Linear 
Model with Non¬ 
stationary Noise 


Figure 5.2 Bark 
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Figure 5.3 Straw 























































(c) Second-order Linear 
Model 


(d) Second-order Linear 
Model with Non¬ 
stationary Noise 


Figure 5.5 Wool 























(a) Original Texture (b) First-order Linear 

Model 



(d) Second-order Linear 
Model with Non¬ 
stationary Noise 


(c) Second-order Linear 
Model 


Figure 5.7 Sand 
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(c) Second-order Linear 
Model 


(d) Second-order Linear 
Model with Non¬ 
stationary Noise 


Figure 5.11 Bubbles 


125 









simulation of sand when the texture is examined in detail 


on a high-resolution display device, 


5.6 Rotation and Magnification 


The effects of image texture rotation and scale 
changes on the covariance function of a texture can be 
determined by expressing each as a linear transformation 
of the linear model. Using the notation of the linear 
model as expressed by Eq. (5.1) we may define 


M 



k=l 


(5.32) 


Using this notation, the maximum likelihood estimators for 
the mean and covariance matrix of N-variate normal 
distribution, N()£:y, T) 


N(*sU,r) = (2^)' N/2 |rr 1/2 exp[-l/2(£-p) T r 1 (^-M)1 (5.33) 


are 


«*-* 


(5.341 


and 


M 


■ sE%v <*- ?i 


(5.35) 


k=l 


Relating this to the facts discussed in sections 3.4, 3.5 
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and ^ , we see that the linear model parameter estimates 
are obtained from the maximum likelihood estimators for 
the mean and covariance matrix. 


If we then define a linear transformation H such that 


U) 


k 



(5.36) 


Then 


M 


to 


= 1 y 


HX, = Hy 


(5.37) 


k=l 


and 


M 


r . = =; 


- k L 


M ^ (HX k ' HX k> (HX k - HX > 
k=l 

- s«IX - v i *£ hT - * ThT > 


k=l 

M 


(5.38) 


- h s XX - x k> <x k - x,ThT 


k=l 


A rp 

= ht h 
-► 

X 

Thus if our model is transformed linearly, the estimates 
required for linear model parameter estimation can be 
obtained from original model estimates and the 
transformat ion itself. 


Rotation may be thought of as a linear transformation 
of coordinate systems. Where our discrete image is 
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considered to be rotated about an axis by angle 0 and the 
standard row, column notation is used (I(n^,n 2 ) represents 
the pixel value for image I at row n^ and column ^) 


I rotated (n l ,n 2 ) 

I . . ,(n_sin6+n,cos0,n,sin6+n„cos0) . 

original 2 1 1 2 


(5.39) 


Usually the row, column addresses in the original image 
are fractional. Therefore, these pixel values must be 
specially defined. A widely accepted practice is to 
estimate each value as a function of pixels surrounding 
it. The most likely candidates are nearest-neighbor and 
bilinear interpolation[42]. It will be shown that in 
either case, the new value may be expressed as a linear 
combination of values in the original image. 


A very similar result may be derived in scale changes 
of a texture. Here 


I i . (n. 
scaled 1 


n 2 ) - 


I . . . (n. ■ 

origin ?l 1 


a,n_*a) 


(5.40) 


The origin of the coordinate system defines the center of 
the image magnification or reduction. In rotation, the 
origin of the coordinate system defines the axis of 
rotation. In both image magnification and rotation, the 
row and column addresses of the pixel in the rotated image 
may be fractional. Again, the value may be expressed as a 
linear combination of values in the original image. 
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Consider the example shown in Fig. 5.12. With angle 
of rotation 0, origin and magnification of 4/3, using 
bilinear interpolation (between rows the columns) we have 


r 


W 1 = [0.65 V 1 +0.35 V 4 ]-0.375+[0.65 V 2 +0.35 V].0.625 

W_ = [0.375 V_+0.625 V c ]•0.35+[0.375 V-+0.625 V,1-0.65 
z zo ,3 b 

W = V 
3 5 

(5.41) 

W 4 = [0.625 V 4 +0.375 V y ]•0.65+[0.625 V 5 +0.375 V g ]*0.35 
W 5 = [0.35 V 5 + 0.65 V g ] -0.625+ [0.35 Vg + 0.65 Vj‘0.375 . 


This implies that 


(5.42) 


0.24375 

0 . 4062b 

0.0 

0 . 1312b 

0 . 2187b 

0.0 

0.0 

0.0 

0.0 

0.0 

0.13125 

0 . 2437b 

0.0 

U . 2187b 

0 . 4062b 

0.0 

0.0 

0 . 0 

0.0 

0.0 

0.0 

0.0 

1.0 

0.0 

0.0 

0.0 

o.o 

0.0 

0.0 

0.0 

0 . 4062b 

0 . 2167b 

0.0 

0 . 2437b 

0 . 1312b 

0.0 

0.0 

0.0 

0.0 

0.0 

0.21875 

0 . 1312b 

0.0 

0 . 4062b 

0 . 2437b 


The accuracy of the estimate depends on the ability oi the 
interpolating function to accurately estimate the value of 
off-grid samples in a discrete image. The covariance 
function itself is being interpolated in this method. 
Caution should be exercised when using a nearest-neighbor 
approach as the transfo r.nat ion of a non-singular 
covariance matrix can be singular if the rows of H are not 
independent vectors. 


By rotating covariance matrices, we are able lo 
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produce models which can generate textures at any angle 
and magnification from one matrix. This could also be 
useful when trying to identify textures of different 
orientation and scale based on covariance statistics. 

5.7 Second-Order Linear Model 


When we 

say that a model is a 

linear or nonlinear. 

we 

are referring to linearity or 

nonlinearity 

i n 

the 

pa rameters. 

The value of the 

highest power 

of 

an 

independent 

variable in the model 

is called the 

order 

o f 

the model. 

For example. 





Y = 6l V 1+ S n B 8 + $„+ 

e 

(5. 

43) 


is a second-order linear model. A general second-order 
linear model with two independent variables may be written 
as 


y =e 1 v 1 +e 2 v 2 +e 11 v^+8 22 v2+e 12 v 1 v 2 +6 0 ♦ E . 


(5.44) 


A full second-order model with N independent variables 
will employ (N*+3N)/2 terms in addition to the 8 

0 

(constant) and e (error) terms. This general second-order 
linear model may be written as 


Vl " 6 iV 6 2 V 2 + '-- +8 nV 6 0 +B 11 V 1 + 6 12 V iV-"* 6 NN V N + E 


N N N 

= 6.V.+ V'' y'' 8. .V.V.+ 8„ + 

/ ^ l x L-t ij l D 0 

i=l j=i 


(5.45) 


i=l 
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Second-order models have been particularly useful in 
studies where surfaces must be approximated by polynomials 
of low order. In all cases, a second-order model will 
"fit" given data as well as or better than a first-order 
model that is a subset of second-order models. This does 
not imply that the second-order model will be more correct 
however, as the process which we are attempting to model 
may be in fact a linear first-order process or some other 
type. 


The use of a second-order model to approximate the 
surface of the general stochastic model could have many 
advantages over a first-order model. An example of 
fitting such a model in one dimension to a given set of 
data is shown in Fig. 5.13. 

Still the linear first-order model may provide a good 


fit to the 

data 

and 

the 

mag nitud e 

of the unexplained 

variance in 

the 

data 

may 

be large 

enough that the 


improvement due to the addition of second-order terms to 
the model may be barely noticeable. In two dimensions, 
the fitting problem is one utilizing a quadric surface 
such as a elliptic paraboloid or hyperbolic paraboloid 
versus a plane to fit a given set of data. Again, the fit 
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may or may not be markedly better. Adding second-order 
terms to a model will always produce a fit as good as or 
better better than a first-order model but the number of 
computations required to compute the coefficients and fit 
the model are much greater. 

It is also important to note that the covariances of 

the are required in order to obtain least-square 

estimates of the parameters 3^ in the first-order model 

[20]. Covariance is essentially a second-order statistic. 

Therefore, estimating the parameters of a second-order 

model will require the use of fourth-order statistics. 

Specifically the correlation of terms V. V. and V. is 

needed. This may cause serious problems as many cases the 

variables in a second-order model will be highly 

2 

intercorrelated. For example, the terms V , and 
(if V ^ is highly related to V-^) may be strongly 
correlated. This situation, often referred to as 
multi-col1inearity, may cause problems during the 
inversion or decomposition of the estimated correlation 
matrix, a necessary step in model parameter estimation. 
For this reason, care should be exercised during the 
analysis of second-order models. 

Inside a circular radius of 14 pixels from V N+ -^ there 
are 307 pixels. To search all possible cross products in 
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this region to find the most significant would require 
over 47,000 cross products to be examined. Computation of 
a covariance matrix containing all of these terms is 
impossible (in practice). In our study we were limited to 
investigate only 820 possible cross products for entry 
into the generation model. As most of the variance was 


ex 

plained 

by t 

he linear terms of 

the model, most of 


the 

cross pro 

due ts 

were 

insig nificant 

from a statistical 

po 

int 

o f 

view. 

This 

se 1 ec 

tion procedur 

e is detailed in [20 

] 

and 

in 

Chapt 

er 3 . 

Tho 

se that were 

significant were en 

te 

red 

into the 

model 

and 

a new texture was generated 

us 

ing 

Eq 

. (5.45 

) wi 

th sta 

tionary Gauss 

ian noise and having 

zero 

mean and 

f ixed 

var ia 

nee [56]. 





The 

resu 

Its o 

f texture 

simulations using 


the 

se 

cond-or 

der 

1 inear 

model are 

shown in Figs. 5.1 ( 

c) 

to 

5. 

11 (c) . 

On 

some 

of these t 

extures only a s 

li 

ght 


improvement from the addition of second-order terms may be 
seen. In most cases, no change can be observed even when 
the results are displayed on a high-resolution display 
device. The lack of improvement could be due to the small 
number of cross-terms examined; however we feel that this 
number is sufficiently large to show any considerable 
improvement due to the addition of second-order terms to 
the linear model. 
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5.8 Textures with Non-Stationary Noise 


Applying a texture generation model to the original 
parent texture image data used to estimate its parameters 
gives a residual error image. When applying the linear 
model to a two-dimensional texture, a two-dimensional 
image containing the pixel differences or residuals 
V„-V„., is found. Here V,,,. is the prediction of the 
next pixel in the sequence as a linear function of the 
pixel around it according to the model without any noise 
added. Naturally, we would expect these errors to be 
small as merely subtracting one pixel from its 
nearest-neighbor would yield a small value in most 
natural, low-noise images. Such an image of residuals was 
generated for the sand and linearly rescaled to show the 
detail present in the image (Fig. 5.14). Definite 
patterns are seen to exist in this image and thus a 
violation of the independent assumption is indicated. 
Ideally, this residual image would be uncorrelated noise. 


A histogram showing the number of V N+ ^ occurring at 
each pixel value is shown in Fig. 5.15. A plot indicating 

shown 

in Fij. 5.16. As would be expected, residuals where the 


the mean of the residuals V„. . -V„. . versus V.,,. i 

N+l N+l N+l 


V N + 1 * s l e53S than 0 will have a mean less than zero and 


those residuals where the V. 


N+l 


is greater than 255 will be 
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likewise positive. Figure 5.17 shows a similar plot of 

A 

the standard deviation of the residuals versus V 

N+l 

These three figures seem to indicate that the distribution 

/v 

of the error in the model is related to the value V N+ ^. 
Therefore the assumption of constant error variance is 


questionable, 


It may be reasonable to drive the 


generation process with noise which does not have 
stationary mean or variance. The effect of such a change 
in the generation process was studied. Figure 5.18 shows 

A 

the distribution of error V N + 1- V N+1 hi stogram of the 
residual image) which appears to be aproximately normal. 


The distribution of this error and the relationships 
between the predicted and actual pixel values was utilized 
to generate textures using non-stationary noise. The 
procedure begins by generating a pixel V N+1 according to 
Eq. (5.45) excluding the error term. With this predicted 
value a random error value e is chosen to be added to 

A 

V N+1* This error value e is chosen from the distribution 

/v 

of error as a function of V and can have any arbitrary 
distribution. The next pixel will than be computed in a 
similar manner. Results of texture synthesis formed using 
this model are shown in Fig. 5.1(d) through Fig. 5.11(d). 


The arbitrary distribution of error as a function of 
V M4 . 7 is calculated by applying the calculated linear model 













Figure 5.17 Standard 

deviation of 

v N+r v N+i vs - 


N+l 


Figure 5.18 Histogram of 

residual image. 
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to the original parent texture and computing a histogram 

/\ 

of errors as a function of , . 

N+l 

In most cases, considerable improvement is seen when 

these .simulations are critically observed on a 

high-resolution display device and compared witli the 

stationary model results. Of course, the information 

required to generate them is considerably greater also. 

The distribution of errors as a function of V must be 

N+l 

condensed and coded to some degree to minimize storage 

A 

requirements. For a 256-grey-level image V N+ ^ usually 
ranges from -50 to 305 and the errors, ^n+i -v n+ 1' ^ rom 
-255 to +255. These ranges were determined 

experimentally. This would yield quite a large amount of 
data if fully stored. By storing a small number (under 
100), typical errors for each range (and not each single 

A 

value) of V the number of data values we are required 
N+l 

to store can possibly be reduced to under 1000. Therefore 
it is believed that this approach of using non-stationary, 
non-Gaussian noise to generate textures may be quite 
acceptable even with severe storage limitations. 

5.9 Conclusion 

The results in this chapter indicate that many 
natural textures are well simulated using a large 


autoregressive model. Adding second-order terms to the 








model improves the results slightly but the resulting 
increase in computational complexity makes this 
second-order simulation method difficult to implement. 
Using non-stationary noise in the generation process 
improves the simulations considerably when the textures 
are viewed critically. The subsequent increase in storage 
and computation required by this addition is small and 
therefore this model modification should be considered in 
most applications. 
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CHAPTER 6 


MULTIPLE MODEL TEXTURE GENERATION 


G.l Introduction 


In this 

chapter, 

we will present three 

methods 

o f 

generating 

textures 

using multiple texture 

models. 

The 

first method 

introduces a set of generation 

kernels 

that 


is used to synthesize texture pixels in a multiple-pass 
manner. Associated with each of these kernels is a unique 
model. This method could be useful in generating textures 
which have very coarse structure. The second method uses 
a piecewise-1inear method of fitting the model to parent 
texture data. The model chosen during the generation 
process is allowed to depend on the pixel values in the 
kernel. Although the fit of the model is better, the 
synthesis results show little improvement over the 
single-model approach with a linear model. The third 
method of texture generation presented in this chapter 
uses a additional image to determine the model number to 
be used during the generation process. This composite 
generation method could be useful when a texture is 
actually composed of a set of subtextures as it allows a 






unique model for each of these subtextures to be used in 
the generation process. 


6.2 Skip-Generate Method 

Simulating textures which have a fine structure is 
usually a much easier process than simulating textures 
with coarse structure. This occurs because the linear 
model contains fewer terms if the texture pixels become 
uncorrelated over a small distance. For the same texture 
at a greater magnification, the pixels become highly 
correlated and the linear model will be forced to contain 
more terms. As the texture becomes more coarse, more 
time-consuming statistical measurements must be taken on 
the parent texture over larger windows. Motivated by 
these problems, the texture generation algorithms in this 
section have been developed. 

In the texture work so far, pixel V was generated 
based on pixels above or to the left of it (see 
Fig. 3.1(b)). As discussed in Chapter 3, the kernel does 
not have to be contiguous. This kernel shape is chosen to 
insure that the image space of our synthesized texture was 
filled during the generation process. However, generating 
pixels along a row, row by row is not the only way of 
filling an image space. 
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Consider the non-contiguous kernel mask in Fig. 6.1. 
If the spacing between the pixels in this mask is 8, using 
the linear model in Eq. (5.4) to generate the right-most 
pixel in the bottom row, we can generate every 8th pixel 
along every 8th row. At each step the next pixel is 
generated based on the previously-generated pixels around 
it (ignoring boundary conditions). After generating an 
image with this type of spacing, the pixels midway between 
the previously-generated pixels on each row may be 
generated using the mask in Fig. 6.2. In this mask, the 
pixel with the "x" in it denotes the next pixel, V ,, to 
be generated according to Eg. (5.4). Naturally, the 
linear model used in this step will have different 
coefficients than the previous one. It is also 
interesting to note that new pixels depend not only on 
previously generated pixels above them (as with the mask 
in Fig. 3.1(b)) but depend also on the pixels below them. 
Still, ignoring boundary conditions, each pixel depends 
only on previously generated pixels. At the next step a 
mask similar to that in Fig. 6.3 can be used to fill in 
the pixels midway between the previously-generated pixels 
in each column. Again pixels are allowed to depend on 
pixels around them. 

By repeatedly using the masks in Fig. 6.2 and 
Fig. 6.3 with successively closer and closer pixel 
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Figure 6.1 First-pass Generation Kernel 
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Figure 6.2 Second-pass Generation Kernel 
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spacing, the texture simulation image space is filled. An 
example showing the pixels generated at each successive 
pass is shown in Fig. 6.4. Wore importantly, to determine 
the linear model for each mask, only one covariance matrix 
is required and can contain as many or as few terms as 
desired. The process of collecting statistics for one 
matrix is not beyond the complexity that we would want to 
undertake for the small number of times required by this 
process. Naturally, any other stochastic process may be 
substituted for the linear model. As before, only the 
measurements required to estimate the parameters 
corresponding to each mask need to be taken. This number 
depends on the spacing of the pixels in the first mask, 
which should be a power of two. Other odd-shaped kernels 
and kernels whose spacing is not a power of two could be 
designed to form space-fi11 ing sets. Most would require 
more models to be estimated and would provide little 
additional information. 

Texture simulations using this method are shown in 
Figs. 6.5-6.12. Only a slight improvement is seen in some 
of the texture simulations over the synthesis done by the 
earlier single linear model. Most of these textures are 
apparently wall simulated by a carefully chosen model and 
the results are not critically dependent on the coarseness 
of the textures. 
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Figure 6.3 Third-pass Generation Kernel 
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Figure 6.4 Filled Space of Skip-generate Method 



Figure 6.5 Skip-generate 
Bark 


Figure 6.6 Skip-generate 
Straw 



Figure 6.7 Skip-generate 
Cloth 



Figure 6.8 Skip-generate 
Wool 







Figure 6.9 Skip-generate 
Leather 


Figure 6.10 Skip-generate 
Sand 



Figure 6.11 Skip-generate 
Water 


Figure 6.12 Skip-generate 
Wood 
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A word of caution should be added concerning the 
computations involved in the linear model coefficient 
calculation of this method. During the later stages of 
the skip-generate method, the pixels in the generation 
kernel become highly correlated as the distance between 
them decreases with each pass. This may cause the 
correlation or covariance matrix of the model to be 
ill-conditioned. To avoid numerical problems, the number 
of variables entered into the process, and therefore the 
number of steps involved in the matrix decomposition 
process, should be kept to a minimum in some cases. The 
use of ridge regression techniques [43,44,46] might also 
be considered. 

6.3 Piecewise-Linear Models 

When generating textures using the general linear 

model described by Eq. (5.4) and the generating kernel in 

Fig. 3.1(b) the same model is used regardless of the 

values of the pixels V^,...,V N . developing more than 

one linear model and allowing the choice of the model at 

each pixel generation step in the synthesis process to be 

dependent on some functional value of V,,...,V„, 

1 N 

F(V,...,) a new synthesis model is formed. 

To illustrate this concept consider the data in 
Fig. 6.13(a). If we were to fit one linear model to the 


149 








data in order to predict from it would look like the 
single line running through the data in Fig. 6.13(a). 
This linear model could then be used to predict V 2 based 
on the value of V . But if we allow the choice of our 
linear model to be dependent on the value of , then for 
an incoming value of we choose a model whose domain 
includes V-^ to predict V 2 . For 5 linear models, the 
straight lines are shown in Fig. 6.13(b). The fit to the 
data using multiple linear models will always be as good 
as or better than that of the single linear model. That 
is, the mean square error will generally be reduced using 
multiple models. 

Using multiple linear models for texture synthesis we 

would generate pixels V based on pixels V ,...V in the 

N+l 1 N 

following way. First, we compute a function, F, of the 
V^,...,V pixels which allows us to choose the proper 
linear model. Then using this model with the values 
V^,...,V N we predict V N+ ^ and add noise. This process is 
diagramed in Fig. 6.14. 

Ideally, the function F should be chosen to minimize 
the total mean square error resulting from fitting the 
limited number of models to the sample data. This is very 
difficult to do in practice however as for W larger than 3 
we are fitting multiple hyperplanes to data in an N+l 
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dimensional space. 


One texture synthesis of sand was done using the 
multiple linear model (see Fig. 6.17). In this case eight 
models were used and the model number was chosen by 
examining the pixel immediately to the left of the pixel 
being generated. The range of this pixel, 0 to 266, was 
divided into 8 equal subranges and the model was chosen 
according to the subrange into which the value fell. Only 
a slight improvement over the single linear model 
synthesis (see Chapter 5) is seen even though the same 
kernel shape was used. No other simulations have been run 
using this model as it is felt that little improvement 
will result. Also the linear model containing 
cross-product terms in Chapter 5 probably provides a very 
good fit in most cases and in more dimensions. In one 
dimension the model of Chapter 6 would fit the data in 
Fig. 6.13 with a quadratic curve. 

6.4 Field-Definition Stochastic Model 

Another method of using multiple stochastic models is 
to generate an image of fields defining the model number 
to be used in a second pass. Such an approach would be 
useful in simulating textures which have multiple 
sub-textures within them. A simple analytical example is 
shown in Fig. 6.15. Real world examples might include 





such things as a brick wall where the texture of the 

bricks is different than the texture of the concrete 
separating them. It was felt that this type of approach 
might be useful in the simulation of bark (see 
Fig. 3.6(a)) which has a strong macro structure. A method 
to separate this texture into two fields, which would 
later define the model to be used, was designed. This 

result (Fig. 6.18) was obtained by successively passing 
smart median filters of varying sizes over a binary image 
(which was obtained by thresholding an original continuous 
grey-level image) (see Fig. 5.2(a)). 

The smart median filters replace the center pixel of 

a window with the median only if certain conditions are 

net. The window is passed over the entire image pixel by 

pixel along a row, row by row in a two-dimensional 

convolution manner as in Fig. 6.15. Let I (n ,n ) be the 

112 

input image and I (n ,n ) be the output. Let the pixel 

0 

values be 0 (black) and 1 (white). Let N denote the 

B 

number of black pixels in the window being processed with 
center anc ^ l et be tbe number of white pixels 

in the window. For the binary case,“the smart median 
filter is defined as 
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Figure 6.13 Linear Models 



Figure 6.15 Analytic 

Sub-textures 


Figure 6.16 Two-dimensional 
Convolution 
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> 
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1, 

if 

W 

V N w 

> 

Thresh 


(n^n^) , otherwise 


( 6 . 1 ) 


Equation (6.1) indicates that the center pixel of the 
window is replaced by the median if the percentage of 
black or white pixels is above a specified threshold. The 
term "replaced" is used loosely in this context to 
indicate the visual appearance of replacement when the 
input and output images are superimposed. In the binary 
case, the median of the pixels in a window is equal to the 
value of the most frequently occurring pixel in the 
window. The threshold and window size varied in the 
successive passes over the image. The window sizes and 
thresholds for each of the passes used to obtain Fig. 6.16 
were 1-0.50, 3-0.50, 5-0.50, 7-0.75, 11-0.78, 7-0.50. 
This multipie-pass procedure helped to retain detail and 
eliminate fields too small for useful measurements. 


The field-definition stochastic model was also 
applied to the non-stationary cloth texture (see 
Fig. 5.14(a)). This texture is clearly composed of two 
alternating subtextures. In this case, the texture fields 
can be extracted by hand (see Fig. 6.21). 
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Figure 6.19 Field-definition Figure 6.20 Field-defini- 
for Bark tion Bark 

Simulation 












Once the two-field images for bark and cloth were 
obtained they were processed again to define a total of 
four fields, the two original fields and two border or 
transition fields. The transition fields are required as 
the kernel of points over which relational measurements 
(usually second-order statistics) are taken may not fall 
wholly within one of the two fields. Such border 
measurements will surely increase the inaccuracy of the 
models for the original fields. So these border regions 
are considered to be unique texture fields. They are 
mathematically defined by passing the original linear 
model texture generation kernel obtained in Chapter 5 over 
the texture and defining the relative importance of each 
pixel to be equal to the absolute value of its associated 
6^ coefficient of the linear model (Eq. (5.4)). These 
coefficients are usually largest near the eye of the 
kernel. The "importance" of the eye itself was set to 
1.5 max(B^). If a large percentage (90% for bark, 95% for 
cloth) of the total "importance" fell within one of the 
two original fields, the point would be considered as a 
member of that field. If not, the point would become a 
member of the transition field based on its position with 
respect to the two fields, A and B. If, at the point 
being processed, A is to the left of B then the point is 
in one transition field, if B is to the left of A then the 
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point is in the other transitional field. The right-left 
orientation is used as the texture generation process is 
usually done in a 1eft-to-right manner along each image 
row. 


The results from processing the original field data 
for the textures bark and cloth are shown in Fig. 6.19 and 
Fig. 6.22. The four fields are represented using four 
grey-levels. The black regions on the border may be 
considered to be undefined as the kernel points do not lie 
with the image region in these areas. 

Tnese field definitions were used both in the 
statistics gathering process as well as the texture 
synthesis process using the four calculated models. The 
models were linear (see Eq. (5.4)) and the methods used in 
their estimation were discussed earlier in Chapter 5. The 
synthesis results are shown in Fig. 6.20 and Fig. 6.23. 
Unfortunately, the field structure is not strongly 
apparent. The primary reason is that each model requires 
a number of pixel generations before it reaches a steady 
state and in most cases this is much greater than the size 
of most of the fields in the bark texture. By observing 
the border regions of many synthesis runs, it seems that a 
steady state is reached after 10-20 pixels have been 
previously generated. In the cloth texture synthesis of 
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Fig. 6.23, the randomness of the noise seems to destroy 
the structure of the texture quite frequently. This is 
illustrated in Fig. 6.24 where only one of the four models 
was used to generate the entire texture. Each of the 
models contained the same set of kernel points used in the 
original linear model of Chapter 5. It is possible that 
improvement could be made by choosing a different set of 
points for each field-model. 

In an actual simulation, once a field image is 
obtained a method must be developed to simulate this field 
texture and this field texture will then in turn be used 
to choose the model numbers in the generation of final 
synthesis. Generating textures with only a few grey 
levels can be done using more complete stochastic 
statistics, perhaps N-grams, but the large size of the 
fields may require that a method such as the skip-generate 
method (discussed in section 6.2) be used. For many 
textures, such as cloth, the field generation process 
would be simple. 

More work should be done in the multiple model area 
as the storage requirements for such models is very small 
but can potentially produce improved results. A great 
number of combinations and approaches are possible in this 
area and, in many cases, may be chosen by the application 
or the particular textures being simulated. 
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CHAPTER 7 


8EST-FIT TEXTURE MODEL 


7.1 Int rod uc tion 

A method of generating texture simulations according 
to their Nth order densities was investigated for binary 
textures in Chapter 3. The simulations resulting from 
this Markov process resembled their parental textures 
quite closely in most cases. When applying a similar 
concept to multi-grey level imagery, the limits of 
computer storage are soon reached. To circumvent this 
constraint, a new method of texture synthesis was 
developed and applied to a number of textures. Simulation 
results using this method are given in this chapter. 

7.2 N-grams in Continuous Imagery 

In binary texture generation based on N-grams a 
single functional value P (V ,/V, , . . . , V„ ) was stored for 
each possible pattern (V^ , , . . . , ) where the \A ' s can be 

zero or one. This value, also called a generation 
parameter, represented the conditional probability that 
the next pixel, , in the generation process would be a 








zero-valued, black pixel. The ' s were chosen by a best 

linear model fit detailed in Chapter 3 and therefore the 

kernel of previous pixels (V^,...,V ) is not necessarily 

contiguous (see Figure 7.1). Details concerning the 

estimation of P (Vjj+l /V 1 ' * * ' ' V N ) ^ rorn a parent texture are 

given in Chapter 3. For binary textures, this single 

value is sufficient to define the distribution of . 

N+l 

given V ,...,V . The number of different functions which 
1 N 

N 

must be stored is 2 . In the generation process each 
pixel is generated based on the values of the pixels 

V ,...,V N surrounding it and on a computer-generated 
uniformly-distributed random variable. The texture 
simulations are generated pixel by pixel along a row until 
each row is complete. Pixel generation along the edges of 
an image can be handled in a variety of ways but in 
Chapter 3 pixels in these border regions were assumed to 
be any random value, 0 or 1, if they were outside the 
image boundaries. 

A similar approach could be used to generate 

multi-grey-level textures. For a texture containing g 
N+l 

grey levels, g different functions, P(V N+ ^/V^,...,V N ), 

roust be stored. (Actually only (g-1) g N are required as 
g-1 

P (X/V^, . . . ,V^)-1 for all V^). Storage limitations 

are soon reached. Also estimation of P(V„ ,/VV ) is 

N+11 N 

difficult as multiple occurrences of the pixel pattern 
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V /•••/V may not exist in the parent texture. Therefore 
1 N 

even without storage limitations tne problems of 

estimating P (V /V ,. ..,V ) from a given parent texture, 
N+l 1 N 

which represents the distribution of V given the values 

N+l 

of V ,...,V is complex. 

1 N 

This estimation problem no doubt has a number of ad 

hoc solutions. The problem is basically that for large N 

and/or large g, there may not be a suitable number of 

occurrences of the pattern V^,... f V to adequately 

estimate the distribution P (V /V ,...,V ) given a finite 

N+l 1 N 

sample size. Even though a certain pattern never occurs 

or rarely occurs in our sample parent texture it is not 

implied that such a pattern is impossible and will never 

occur in our simulation synthesis. We might often find 

numerous occurrences of this pattern if our sample size or 

the size of our parent texture was increased, especially 

in noisy and fine-structured textures. But as this very 

large sample may not be present, we must estimate 

P(V /V ,...,V ) for all V ,...,V based on available 
N+l IN IN 

samples. 

One approach would be to use sample patterns which 

closely resemble but which may not be exactly the same as 

each pattern (V ,...,V ). That is in a pictorial sense, 
1 N 

we use patterns of (V_,...,V ) which look "close to" the 

1 N 





pattern for which we are attempting to estimate 

P(V /V f ••-#V^). Therefore samples in our sample parent 

texture may be used to estimate numerous P (V.,,,/V, .. •. ,V„) 

N+l 1 N 

and not just those they fit exactly. The concept of a 

distance function must be used to numerically define 

"close to". Given two patterns, one from our sample 

texture and the other from the conditional probability of 

the kernel we are attempting to estimate, the distance 

measure can be used to determine the value of that sample 

in estimating P(V /V ,...,V ). If the fit between the 

N+l 1 N 

kernel pattern and the pattern in the sample texture is 

good the associated value of V^ +1 in the parent texture 

will be valuable in estimating P (V /V ,...,V ). 

N+l 1 N 

Normally, when N and g are small or when we have many 
samples for any given v g'-*-» v N ' we can use the histogram 
of the associated V N+ j_ to estimate P (V N+ ^/Vj,...,V N ). 
Here the relative number of times a particular value of 
V n+i occurs given a pattern indicates the conditional 
probability we are attempting to estimate. This was 
discussed in section 3.2. tvnere a distance measure is 
used instead, a good fit could be considered to be 
synonymous with high frequency of occurrence of that 
pattern and a poor fit with low frequency of occurrence. 

If such a method of estimating these conditional 
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probabilities is used we are still faced with a huge 
storage problem. For this method to be practical, the 
storage requirement n.ust be reduced. From an information 
standpoint, it is interesting to note that a method of 
estimating N-grams or conditional probabilities 
P(V N+ i/Vi,...,V N ) from a sample parent texture image 
produces g N+ ^ data values from M pixel samples where M is 
the size of the square parent texture image in pixels. 
For large g and N this is a drastic increase in data. But 
the actual information content can really never be greater 
than that content of the sample parent texture image. 
Therefore, this M value represents an upper bound on the 
amount of data we should use to generate a texture 
simulation. Any amount of data exceeding this will 
contain redundant data. 

7.3 The Synthesis Method 

Combining this concept of upper bound with the idea 
of forming a distance measure to compare two texture 
kernel patterns leads to a new texture synthesis method. 
In this method, we generate the next pixel based on the 
pixels in the kernel surrounding it (see Figure 7.1) and 
their comparison to similar kernels in the parent texture. 
This comparison is made by passing the kernel currently 
present in the simulation process over the parent texture 
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and computing the distance function at all possible points 
(see Figure 7.2). Denoting the pixels in the parent 
texture by X , i,j=0,...,/M-l and the pixels in the 

i/j 

kernel V.,...,V by Y. . , we can compute a comparison 
-L N 1 r J 

image 

c , = COMPARISON (X . , . ,Y. .) , n 

a,b i+a,g+b' i,] (7.1) 

for all a,b such that the kernel is within the boundaries 
of the texture. 

One possible comparison function would be 

correlation. Assuming, without loss of generality, that 

our kernel is contiguous as in Figure 7.3 and the elements 

are denoted as Y. ., this function would be defined as 
1 > 3 



The problem with this particular distance measure is quite 
serious. Correlation does not take into account 
differences in over-all mean. For example, the kernels in 
Figure 7.3 are perfectly correlated but their means differ 
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significantly. Differences in first-order statistics 
between these kernel patterns will not be detected by a 
correlation measure and so another comparison function to 
supplement a correlation value would be required. 

A better comparison function would be the mean square 
difference (MSD). This is defined as 

MSD , =) / (X.^ . , —Y. ,) 2 

a,b L—iL—i i+a,g+b (7.3) 

i j 

where i,j must be within the coordinate range of the 
kernel as in Eq. (7.1). This comparison function will 
detect differences in first-order statistics between two 
pixel patterns (such as those in Figure 7.3) as the MSD 
function is a sum of squares of differences. Whereas the 
correlation coefficient of Eq. (7.2) varies between -1.0 

and 1.0 and is largest when the fit is good, tne MSD 
measure is small when the fit is good and it is always 
positive. 

The MSD function weights the comparison of all 
elements in a kernel equally. Having studied many texture 
generation models we immediately recognize that this fit 
is not properly weighted. The few pixels which are 
closest to Y next in proximity are far more important when 
predicting y next t * ian those which are far away. So 




Eq. (7.2) must be modified to give the weighted 
mean-square difference (WMSD) 

WMSD , —Y. . ) 2 . W. . . (7.4) 

a,b i+a, j+b i,j' i,j 

i j 

A possible choice for W is 


W. 


(7.5) 


1 / D 


2 2 
* i-i NEXT^ + ^“^NEXT^ 


where R is the euclidean distance between pixel Y. . and 

1/3 

the kernel eye Y„_„ m and the coordinates of the eye are 

NEXT 1 


given by (i NEXT -3 NEXT > - 


As the first step in comparing a given kernel Y^ ^ to 

all kernels in the parent texture, for each point (a,b) in 

the parent texture, ignoring edges, the WMSD is computed 

resulting in an image of WMSD's. Where the fit between 

the generated kernel Y. . and the image X. .is good, we 

1 / J 1 I D 

would expect WMSD^ ^ to be small. The smallest WMSD 
represents the "best" fit according to our norm. We could 
choose the Y jjEXT assoc i a ted with this best fit at point 
(a,b) to be our next pixel in the generation process, 
however this can cause problems. First of all, the 
generation process would "lock in" on the parent texture 
and the generated texture could very well become just an 
exact copy of the input parent texture. Second, we know 
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p 


ideally that Y has a distribution, not just a mean. 

NEXT 

In the autoregressive model of Chapter 5 we gave Y NEXT a 

distribution by adding random noise to it. Although this 

could be done here, such an approach would fail to use 

additional information contained in the WMSD image. There 

may be a set of points (a,b), all exhibiting a good fit to 

the kernel pattern Y. ,. In fact, the best fit may have a 

i/3 

noisy Y and the other good fits could provide 

NEXT 1 

information to improve the prediction of the Y NEXT in the 
generation process. Using a set of best fits is 
equivalent to increasing our sample size. We look at a 
set of similar patterns to pick our Y 

NEXT 


At 

this point there are 

numerous 

ways 

to proceed. 

Logically those patterns 

with the 

" best 

" fit should 

provide 

better estimators 

for y NEXT 

so 

some kind of 


weighting decision is needed to choose the relative 

importance of the WMSD's found. If we search through the 

WMSD image and find the minimum value, WMSD . , and scale 

min 

all the WMSD's by that we form a new image MAXI 


MAXI 


a, b 


WMSD . 
min 


WMSD 


a,b 


(7.6) 


This image has the value 1.0 at the best fit point and 
values 0 < MAXI < 1.0 elsewhere. 
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Here we can look at the MAXI (a,b) image and study its 

range. I£ 0.16 MAXI £ 1.0 it is implied that the worst 

fit yields a 0.16 value. Somehow that worst fit should be 

translated to imply that the probability of choosing the 

Y NEXT assoc i atec 3 with that point (a,b) W0RST i s nearly 0.0. 

The simplest way of doing that is to take powers of the 

image MAXI (a,b). The maximum remains 1.0 while smaller 

numbers approach 0.0. For example (1.0)^=1.0 but 
10 -8 

(0.16) =1x10 . We do this to obtain an jad hoc estimate 

of P ^ Y NEXT //Y i j)* After experimentally studying the 
values of MAXI (a,b) and its powers, the value of 16 was 
chosen and a new image PDFUNS 

PDFUNS , = (MAXI . ) 16 

a,b a,b (7.7) 

was used to estimate the probability density function 
P (Y xti:iv .„/Y . .). The values in the PDFUNS image are 

Nila I 1 1 J 

generally very small with less than 1% of the image points 
having value greater than 0.1. As a rule of thumb, it can 
be argued that 1% to 0.05% of the 128x128 PDFUNS values 
should be greater than 0.1. A larger percentage would 
increase undesired randomness and noise in the synthesized 
image and a smaller number could cause "lock in" on the 
parent texture. The value 16 was also chosen for 
convenience and computational efficiency as it can be 
computed with only 4 multiplications and minimal data 
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storage. More study could be made concerning the effect 
of the value of this variable on the simulation results 
and other approaches to creating a PDFUNS image from the 
MAXI values could be tried. Studying a histogram of MAXI 
values might also be very informative. 


PDFUNS is then scaled so that 



a b 


PDFUNS (a ,b)=1. 


In this way a pseudo-density function is formed. Finally 
a uniformly distributed random variable, r, [0,1] is 
generated and a point (c,d) is found such that 


c-1 


d-1 


EE' CF “ S a,b* E 

a=l b b=l 

d d 

PDFUNS 


PDFUNS , < r 

c, b 


(7.8) 


EE 

a=l b 


, b + E PDroNS c,b > r • 


b=l 


The Y associated with the kernel shape at (c,d) is 
then used as the next pixel in the generated image. The 
process is continued until a full texture image is 
generated with the kernel window moving one pixel at each 
step, row by row. 


In an indirect way, this is equivalent to generating 
a random variable having any distribution using the 
desired cumulative distribution combined with a uniformly 
distributed random variable (which is easy to generate). 
In other words, uniformly-distributed deviates are 
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transformed to deviates having the desired distribution 
using the inverse cumulative density function [45,48,58]. 
This is frequently done in simulations. 

7.4 Results 

For a kernel containing 55 pixels (see Figure 7.4) 

passed over a 128x128 parent texture approximately 

7.2xl0 6 operations (additions or subtractions) are needed 

to get the WMSD image defined by Eq. (7.4). Another 

2.6xl0 5 are required to find the next pixel according to 

Eq. (7.8). therefore, to generate a 512x512 texture 

12 

requires 1.96x10 (2 trillion) operations. 

Results from texture synthesis done by this method 
are shown in Figure 7.5 through 7.15. The original 
textures are shown in Figs. 5.1(a) through 5.11(a). Each 
of these images is 512x512 pixels. A 128x128 section of 
each original (parent) texture was used for the 
simulation. Bark exhibits very large macro structure and 
this is lost in the simulation. A similar thing happens 
with raffia as the kernel size is smaller than the cell 
size of the original texture but is not as pronounced. 
The top part of the bubbles texture was generated using a 
128x128 portion different than that of the bottom part. 
For this reason the top 20-30% of the texture looks 
different from the rest. The large number of operations 


171 






.1 Two-dimensional 
Synthesis Kernel 


Figure 7.2 Passing Kernel 
Over Parent 
Texture 



51 

52 

53 

56 

57 

60 

63 

66 

*52 

55 



Figure 7.3 Perfectly Correlated Kernels 
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Figure 7.4 Best-fit Model Kernel 


























Figure 7.5 Best-fit Grass Figure 7.6 Best-fit Bark 



Figure 7.7 Best-fit Straw Figure 7.8 Best-fit Cloth 
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Figure 7.9 Best-fit Wool Figure 7.10 Best-fit 

Leather 



Figure 7.11 Best-fit Sand Figure 7.12 Best-fit 

Water 
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Fiqure 7.13 Best-fit Wood Figure 7.14 Best-fit 

Raffia 



Figure 7.15 Best-fit Bubbles 
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makes this process very time consuming even when 


a 


pipelined processor is dedicated to the task. About 5.5 
days of dedicated time on an AP120B were required to 
generate each texture. 

Although this method is of little practical use due 
to the computational complexity of the algorithm a few 
points should be made. With constantly increasing 
computer processing speeds, a simplified version of this 
texture simulation method may be implemented in the near 
future. It is even possible that such computations could 
be performed by an array of microprocessors. In any case 
such brute-force approaches are simple and could be made 
cost-effective in the future. 

The results also indicate visually the amount of 
texture information present in a 55 pixel window (see 
Figure 7.4) because at each pixel generation step, the 
next pixel in the Markov process depends on only the 
pixels in this neighborhood. 

Finally, this approach is admittedly ad hoc . 
Numerous distance measures could replace the one chosen in 
this work and each would give different results that might 
appear better or worse. It is always important that the 
process be random and not merely copy the texture sample. 
If the simulation region is much larger than the parent 
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sample, a deterministic process will quickly generate 
patterns that can easily be seen to repeat. In other 
words, the histogram represented by p ^ V n+ 1^ V 1' * * *' V N^ 
should rarely be a delta function. A reduction in the 
number of computations could be made if the kernel was 
non-contiguous. Also, better results could probably also 
be obtained if the kernel window was larger. The shape, 
contiguity and size of the kernel in this study was chosen 
primarily for computer programming considerations. 

7.5 Conclusions 

The results from this best-fit texture synthesis 
method are very pleasing but the number of computations 
required is large. Other similar algorithms could be 
developed which are simpler and could possibly produce 
even better results. With the decrease in computation 
costs and the increase in processor speeds of future 
computers, such texture synthesis methods could be 
implemented in the future without great cost. 
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CHAPTER 8 


NON-HOMOGENEOUS TEXTURES 

8.1 .Introduction 

In this chapter, methods for removing and introducing 
non-homogeneities in texture images are presented. 
Non-homogeneities in neighborhood mean and standard 
deviation are often removed previous to statistics 
collection over a parent texture image to improve the 
accuracy of parameter estimation. Similar 

non-homogeneities may be added during the texture 
synthesis process by merely reversing the process. In 
this way, synthesized textures which are homogeneous may 
be processed to be non-homogeneous. 

8.2 Removing Non-Homogeneities 

Prior to simulation attempts, the textures in this 
study have been preprocessed by statistical differencing 
[42]. This preprocessing step is described by 

I(n.,n 9 ) = [F(n ,n )-F(n. ,n_) ] -—- + 

2 Ua ( .vn 2 f„ d J (8-1) 

[am d +(l-ra)F(n 1 ,n 2 )] 
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where and represent desired mean and standard 
deviation. F is the input pixel at location (n^,n 2 >» row 
n^ and column n^ in the discrete digital image matrix, and 
I is the output pixel in the statistical differenced 


image. 

F (rif , n 2 ) and 

S ( n i, n 2 

) represent 

the 

mean and 

standard deviation of 

the 

input 

image 

at 

the point 

(n 1 ,n 2 ) 

. The variable 

A i s a 

gain 

factor 

that 

prevents 

overly 

large output 

v a 1 ue s 

when 

a is small. 

and a is a 


proportionality constant controlling the extent to which 
the mean of the output image is homogeneous. 

In our studies, mean and standard deviation factors 
were computed in non-overlapping 16x16 pixel blocks values 
are used to compute the mean and standard deviation at 
each point. In our work a = 0.8, m^ = 128, 35 and 
A = 6. These values tend to induce local homogeneity in 
mean and standard deviation over an image. Large amounts 
of variation, however, will only be reduced and not 
eliminated unless A is very large and cx = 1.0. 

A 

Assuming that o does not approach zero, then another 
form of statistical differencing can be used. This may be 
written in equation form as 
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( 8 . 2 ) 


r 6 °d 

I(n 1 # n 2 ) = [F(n 1 ,n 2 )-F(n 1 ,n 2 )] —-+ 

L 6 (n^,n 2 ) 

+ [am d +(l-a)F(n 1 ,n 2 )] 

Here, both a and 6 are proportionality constants ranging 
from 0.0 to 1.0. Setting a. = 1.0 and 6 = 1.0 returns an 
output image with precise desired mean and standard 
deviation. Setting a = 0.0 and $ = 1.0 causes the 

standard deviation but not the mean of the input image to 
change. Setting a = 1.0 and 6 =0.0 causes the mean but 
not the standard deviation of the input image to be 
modified. Setting a = 0.0 and 6 = 0.0 produces no change. 
Examples of statistical differencing are shown in 
Figs. 8.1 through 8.4. The cork texture of Fig. 8.1 is 
non-stationary in mean due to shading differences, 
primarily at the right edge. Figure 8.2 shows the image 
resulting from processing Fig. 8.1 using the statistical 
differencing algorithm. The non-homogeneity of mean is 
removed and the contrast is slightly increased. An 
original brick texture image shown in Fig. 8.3 has very 
low contrast. After statistical differencing, local 
contrast is much improved and the texture is more apparent 
(see Fig. 8.4). Thus the statistical differencing 
algorithm is quite useful in eliminating 

non-homogeneities. 


(1 


- 6 ) 
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Figure 8.1 Before Statistical Figure 8.2 After Statisti- 
Differencing cal Differ¬ 

encing 



Figure 8.3 Before Statistical 
Differencing 


Figure 8.4 After Statisti¬ 
cal Differ¬ 
encing 


181 









8.3 Introducing Non-Homogeneities 


The inverse operation of statistical differencing can 
be called local moment modification. Solving for F in 
terms of I using Eq. (8.1) we find the formula for local 
moment modification as 


Ac? (n. ,n )+o F(n 1 ,n, ) )Aa 

F( ni ,n 2 ) = -±— -- I(n lf n 2 )+ — -— 


Aa 


- (am d +(l-a)F(n 1 ,n 2 )] 


Ao (n^,n 2 )+a d 


(8.3) 


Using Eq. (8.3) we can introduce non-hornogenei ti es into a 
simulated texture by generating an image F(n^,n 2 ) and 

-A 

o(n lf n 2 ) and defining A, a, n? d , and . 

Again, if we assume that a(rii ,n 2^ does not approach 
zero, then Eq. (8.2) can be inverted to form another local 
modification formula 


F(n d ,n 2 ) = F(n lf n 2 ) + 


o(n 1 ,n 2 )[I(n 1 ,n 2 )-am -(l-a)F(n 1 ,n 2 )] 


(8.4) 


a (n 1 ,n 2 )(l-<5)+6c? d 

In the process of local moment modification, it is best to 
set m^ and to be equal to the mean and standard 
deviation of the homogeneous texture I(n^,n 2 ). Then, the 
mean and standard deviation of the output image (F(n^,n 2 ) 
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1 

1 

! 

will be defined by the images F(n] L ,n2) and a(n^,n2). 

These images may be generated randomly. 

An image, before and after local moment modification, 
is shown in Fig. 8 . 5 . Here F(n^,n2> was assumed to be 

A 

ramp-like and o(n^,n2) was constant. Many other complex 

_. /v 

and random F(n^,n2) and afn^nj) images could be used to 
create different effects and simulate phenomenon such as 
non-homogeneous lighting effects. 
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CHAPTER 9 


TEXTURE IDENTIFICATION AND SEGMENTATION 

9.1 Introduction 

In this chapter we examine various approaches to 
texture segmentation and identification using the linear 
model developed earlier. These methods employ covariance 
measures of a texture over windows of the region being 
segmented or identified. These same covariance measures 
were used to estimate the linear model parameters which 
were examined in Chapter 3 and Chapter 5. In Chapter 5, 
the information content of these measures was shown by 
performing simulations of various textures and therefore, 
based on these results, we might conclude that these 
second-order statistics would be useful in the 
segmentation and identification of textures. Pictorial 
segmentation results are given in this chapter. 

It is generally agreed that a great portion of 
texture information is contained in the second-order 
statistics of a texture. There are notable exceptions to 
this rule as was shown in Chapter 2, however for most 
natural textures, second-order statistics have proved to 
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be sufficient to adequately discriminate textures in most 
applications [49,51]. In a practical sense, we usually 
are also prevented from gathering and analyzing 
higher-order statistics because of computational 
limitations. In fact, we can easily be overwhelmed by 
masses of data arising from second-order statistics. 

The most general second-order statistics are the 

complete second-order joint densities, or 2-grains, which 

were discussed in Chapter 2 for the one-dimensional 

texture synthesis case. These measures may be estimated 

by joint gray scale histograms taken over a parent 

texture. For a g grey-level image, each histogram 

2 

requires g storage locations. But, as was pointed out in 

Chapter 2, there is one such histogram for each vector 

distance or spacing between a pair of pixels. To compute 

these histograms over all spacings, 2,...,n , in two 

2 2 

dimensions would result in (n r -l) g entries. Even for 
reasonable g and n , this could easily create a data 
expansion containing unnecessary information rather than a 
data reduction yielding measurements with discrimination 
value. 

For these reasons, texture image data is often 

quantized (to reduce g) and second-order measurements 
resulting in joint grey-scale histograms are made over a 
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small number of pixel spacinys (to reduce n r ) . To 
decrease the size of the feature space further, various 
functions of joint yrey-scale histograms are calculated 
and these values are primarily used to identify a texture. 
For a single nistogram, as many as 25 to 30 functions have 
been proposed [7], In spite of the large dimensionality 
of the feature space and the problems with quantizing low 
contrast textures, this family of texture features is used 
frequently to successfully classify textures [7,52]. 

Many other identification and classification schemes 
exist [49] as the discrimination of textures represents 
the most important application of texture analysis. 

In this chapter, we reduce the information contained 
in the joint grey-scale histogram to one single number, 
the correlation coefficient for that particular pixel 
spacing. It is expected that this large reduction will 


cause a decrease in discrimination power 

as the 

size 

and 

information 

content 

of the feature 

space 

has 

been 

significantly 

reduced. 

The purpose of 

this exercise 

is 


not to develop a new, more powerful texture identifier but 
merely to access the information content of the 
correlation coefficient values when applied to the problem 
of texture discrimination. It is already apparent from 
the simulation results presented in earlier chapters that 
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a great deal of texture information may be obtained by 
proper use of these correlation coefficients. 

Correlation measurements have been applied in various 
ways to the problem of texture identification previous to 
this study [2,11,12,14,49,53]. It has generally been 
concluded that they are of lesser value than Haralick's 
family of functions on values of the joint grey-scale 
histogram [51] when applied to the discrimination problem. 
In the remaining sections of this chapter we will develop 
two new discrimination methods utilizing correlation 
values. One is based on the statistical test for equality 
of covariance matrices. The other utilizes multiple 
statistical tests for the equality of individual 
correlation coefficients. Both show good discrimination 
power but neither exceeds the quality of Haralick's 
measures. 
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rk 

, measuremen 

ts were 

taken 

ov 

er a 


square, W BIG pixels in length, and a center square, Wgjk^^LL 
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pixels in length, was classified from these measurements 
as in Fig. 9.1. 

Second-order statistics must be measured over a 
variety of vector distances (pixel-pair spacings) to be 
useful in texture discrimination. In our case, these 
second-order statistics are correlation values. There is 
only one correlation value for a particular spacing. 

There are many approaches to estimating a correlation 
value over a window for a particular pixel-pair spacing. 
One involves passing a kernel of pixels over the window 
and taking a sample at all points where the kernel is 
completely contained within the window boundaries. Such 
measurements would result in a covariance matrix for the 


kernel ov 

er that window. 

This 

matrix 

can 

be 

used 

to 

identify 

a texture, as 

will 

be sho 

wn 

later 

in 

this 

section. 

Another approach 

to 

estimati 

ng 

a coi 

'rela 

tion 


value for a particular pixel-pair spacing would involve 
measurements over all possible samples within the window 
of that spacing. This is equivalent to passing a kernel 
containing two points over the window for each pixel-pair 
spacing. The result of this approach is a correlation 
value for each spacing which must be examined by itself 
and may not be used to form a correlation matrix as was 
discussed in section 3.7. A method for identifying 
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texture using these individual correlation values will be 
examined in section 9.3. 

Once we have obtained a covariance matrix by passing 
a kernel of points over a window, that matrix may be 
compared statistically with covariance matrices from known 
textures. Such a statistical test has been derived for 
testing the equality of covariance matrices. The maximum 
likelihood ratio approach used to derive this test makes 
the standard assumptions of multi-dimensional normality 
and independence of samples. Both were used in our 
texture simulation work. Details concerning the test and 
its derivation are given in [47,54]. The statistical test 
for two covariance matrices is given by 

°1 “ 2-3026 dD S X 2 n(n+1)/2 19.1) 

where 

D = (M 1 +M 2 -2)log 10 |c|-(M 1 -l)log 10 |C 1 |-(M 2 -l)log 10 |C 2 |. 

C-^ and C 2 are the estimated covariance matrices for each 
g roup, 

C = [(M 1 -l)C 1 +(M 2 -l)C 2 l/(M 1 +M 2 -2) 

and 

d - 1 - [(M^ry + (M^TI - Tm^ 2)][ ,2n2 - 3n - 1 » /6in+1 >_- 

l\| is the size of the covariance matrices being tested 
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which is equal to the number of pixels in the kernel. 
and M 2 are the sample sizes used to estimate and C 2 . 

|C| denotes the determinant of C. The test statistic, , 
has an approximate chi-square distribution with N(N+l)/2 
degrees of freedom and approaches 0 as approaches C 2 . 

Having derived the test for equality of two 
covariance matrices we must determine the contents of 
these matrices. As the number of points in the kernel 
increases, the size of the covariance matrices increases 
leading to more difficult and time-consuming determinant 
calculation required in Eq. (9.1). Also, as the spacing 
of these pixels increases, the number of samples in a 
window decreases. Finally, as more and more points are 
included in the kernel, the amount of redundant and 
overlapping information in the covariance matrix increases 
due to redundant pixel-pair spacing. This was discussed 
in section 3.7. 

To eliminate this redundancy, we will consider the 
problem briefly in one-dimension. Perfect, non-redundant 
pixel spacing is possible only for patterns with maximum 
range of 1, 3 or 6 pixels given by the corresponding 
patterns XX, XX-X and XX—X-X as shown in Table 9.1. For 
these three particular ranges, the patterns shown are 
constructed so that no two pairs of X's are separated by 
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redundancies 


contain 


if all possible distances and 


orientations less than or equal to the maximum range are 


to be spanned. Patterns obtained using this vector 


product approach are heavily weighted in the horozontal 
and vertical directions. 


A two-dimensional kernel containing 16 pixels which 

spans a two-dimensional range of 6 is shown in Fig. 9.2. 

Passing this kernel over a window yields a covariance 

matrix of dimension 16, which is not an unreasonable size 

for computation purposes. The number of data samples over 

2 

a window of size W is (W -6) . Even for small 

BIG BIG 

windows of size 16, a reasonable sample size of 100 may be 
obtained. 


Once the kernel and procedure are determined, we 
proceed with the process of segmenting the textures to 
test the identification procedure. There are many 
approaches to this problem which are usually defined by 
the particular case of interest. We may or may not have 
prototype and parent texture data. We can segment, 
cluster-analyze or identify. Notice that the generality 
of the test leaves a wide variety of options open. We can 
compare matrices in one area of an image with those in 
other areas and our test statistic values will be the 
distance measures defining the closeness of the textured 
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Figure 9.1 Segmentation Window 



igure 9.2 Two-dimensional Spanning 
Kernel 









Table 9.1 MINIMAL SPANNING SETS 


RANGE 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 
21 
22 
33 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 

44 

45 

46 

47 

48 

49 

50 

51 

52 

53 

54 

55 

56 


IPTS 

2 XX 

3 XXX 

3 XX-X 

4 XXX-X 

4 XXX—X 

4 XX—X-X 

5 XXXX-X 

5 XXX—X—X 

5 XXX-X—X 

6 XXXX—X-X 

6 XXXX-X-X 

6 XXXX-X-X 

6 XXX-X-X—X 

7 XXXXX-X-X 

7 XXXXX-X-X 

7 XXXX-X-X-X 

7 XXXX-X-X-X 

8 XXXXXX-X-X 

8 XXXXX-X-X-X 

8 XXXXX-X-X-X 

8 XXXXX-X-X-X 

8 XXXX-X-X-X X 

8 XXX-X-X—X—X-X 

9 XXXXXX-X-X-X 

9 XXXXXX-X--X-X 

9 XXXXX-X-X-X-X 

9 XXXXX-X-X-X-X 

9 XXX-X-X—X—X—XX 

9 XXX-X-X—X—X--X-X 

10 XXXXXX-X-X-X-X 

10 XXXXXX-X-X-X-X 

10 XXXXXX-X-X-X-X 

10 XXXXX-X-X-X-X-X 

10 XXXX-X-X-X-X—X—X 

10 XXX-X-X— X— X— X— X-X 

10 XX-X—X-X-X-X-X-XX 

11 XXXXXXX-X-X-X-X 

11 XXXXXX-X-X-X-X-X 

11 XXXXXX-X-X-X-X-X 

11 XXXXX-X-X-X-X-X-X 

11 XXXX-X-X-X-X-X-X-X 

11 XXXX-X-X-X-X-X—X—X 

11 XX-X—X-X-X-X-X-X XX 

12 XXXXXXX-X-X-X-X-X 

12 XXXXXXX-X-X-X-X-X 

12 XXXXXX-X-X-X-X-X-X 

12 XXXXXX-X-X-X-X-X-X 

12 XXXX-X-X—X-X-X X XX 

12 XXXX-X-XX-X-X-X-X—X 

12 XXXX-X-X-X-X-X-X—X—X 

13 XXXXXXXX-X-X-X-X-X 

13 XXXXXXX-X-X-X-X-X-X 

13 XXXXXXX-X-X-X-X-X-X 

13 XXXXXXX-X-X-X-X-X-X 

13 XXXXXX-X-X-X-X-X-X-X 
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regioi We could compare a matrix with a number of known 
prototypes and classify the unknown region according to 
the test statistic for each comparison. Finally, we may 
merely wish to locate a particular texture in an image and 
thus we would compare matrices with only one prototype. 

This later approach was selected in this chapter as 
the visual results are simple to display and analyze. We 
will compare the matrix measured over a large window with 
one obtained over the complete parent texture prototype 
and assign the pixels in the small window accordingly (see 
Fig. 9.1). Figure 9.3 shows the test composite image 
used. We will attempt to identify the texture sand in 
that region. The texture sand was chosen as the 

simulation results of this texture were in some ways very 
poor when examined critically on a high resolution device. 
Therefore, it represents a worst-case example. 

Figure 9.4 shows the segmentation results when 

W = W = 16. The pixel brightness are linearly 

mapped test statistic values, which are supposed to be 

chi-square-distributed. Improvement is made when 

W =32 and W =16 as is seen in Fig.9.5. 

BIG SMALL 

Figure 9.6 shows much improved results when W = 48 and 

* BIG 

w„ = 16. The difference is due to the availability of 

SMALL 1 

more texture information in the larger window in the form 








Figure 9.3 Input Composite 
Image 


Figure 9.4 Segmentation Using 
Covariance Matrix, 



Figure 9.5 Segmentation Using Figure 9.6 Segmentation Using 
Covariance Matrix, Covariance Matrix, 


W 


BIG 


32 


W BIG = 48 
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of increased sample size (100 to 1764). 

As a final note, it should be understood that this 
chosen texture kernel (see Figure 9.2) totally ignores 
information obtained earlier concerning the 
interdependence of pixels in the generating kernel as 
expressed in the linear model. It is very possible that 
better segmentation results could be obtained by using the 
kernel shape which best fits each texture as this 
particular pattern of points was found to be most 
significant. The linear model used for simulation of each 
texture might even be used itself as the distributions of 
the parameter estimates are known if certain assumptions 
are made [24,39]. But neither of these approaches have 
been tried even though improvements are expected. 

9.3 Segmentation Using Individual Correlation 
Coefficients 

Individual correlation estimates for particular 
pixel-pair spacings may be made by taking measurements 
over all possible samples within the window containing 
that spacing. This approach utilizes more information 
than the fixed kernel method as it includes measurements 
near all edges of the window (thus, the sample size is 
increased). The result is a single coefficient for each 
pixel-pair spacing. If the spacing between each pair is 
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unique then the set of coefficients will be non-redundant. 
The covariance matrices used to segment textures in the 
last section may contain redundant information if the 
pixel-pair spacings in the kernel are repeated. 

A statistical test for comparison of two correlation 
coefficients has been derived for the bivariate normal 
distribution and is 



where 

U 2 = (Z 1 -Z) 2 (M 1 -3)+(Z 2 -Z) 2 (M 2 -3) 

and 

Z^ = arctanh (r^) 

Z 2 = arctanh (r 2 ) 

Z 3 = [(M 1 -3)Z 1 +(M 2 -3)Z 2 l/(M 1 +M 2 -6) 
r^ and r 2 are the computed correlation coefficients from 

the two groups. This test is derived using the fact that 

is approximately distributed N(Z-^: arctanh (P-^) , 

-1 

(Mj_-3) ) where is the true correlation coefficient of 

the population (see section 3.5). For additional details 
see [39]. Again, to derive the tests assumptions of 
normality and independent samples are made. 

As in the previous section, the test statistics may 
be used to segment or identify textures. However, 
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difficulty arises as our test must be performed for each 


measured correlation value. Thus a number of U 2 values 
are obtained if multiple pixel-pair spacings are used. At 
this point these values can be combined in numerous ways. 
As each U 2 value is chi-square distributed we can compute 

the probability, p n , that a random variable, which has 

u 2 

that chi-square distribution, is greater than or equal to 
M 2 as shown in Fig. 9.7. For U 2 = 0 (r-^ = r 2 ) this 
probability is 1.0. A product of the probabilities 
associated with each U 2 was used to indicate the overall 
fit of one set of correlation coefficients to another set. 

Finally, it should be noted that tests for equality 
of correlation coefficients will not detect differences in 
mean and variance over the window. The correlation 
coefficient specifically removes these effects. As a 
result, it may be advantageous to include tests for 
equality of means and variances into the comparison 
process. For equality of means the test statistic is 


» 1~ V 2 


/(M 1 -i ) ^(M 2 -i)lf JiTT 

’ M, + M- - 2 * M 1 


= t M 1 +M 2 -2 


(9.3) 


and o 


where y, and y 9 are the calculated mean and a _ 

±4 12 

are the calculated variances for the two groups being 
compared. This test statistic is a value of a random 
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variable having the t-distribution with M^+M 2~2 degrees of 
freedom. For the equality of variances the test statistic 
is 


t-, = 


1 

Sf 1 Mi-l.Hj-i 


(9.4) 


This test statistic is a value of random variable having 

an F-distribution with M^-l and N^-l degrees of freedom. 

The derivation of both of these tests requires the 

assumption of normality for the two groups and an 

independent sample set. The known distributions of the 

test statistics t and f^ was used to obtain the 

probabilities, p and p_ that a random variable having 

1 

that distribution would be greater than or equal to the 
calculated values of 11 ^ | and MAXCf^l/f^) or less than 
— 11 ^| and MIN (f ^ , 1/f ^) (see Figs. 9.8 and 9.9). 


Having obtained the set of probabilities, p , 

U 2 


and 


p and p f corresponding to the set of tests for the 
fc l 1 

equality of correlation coefficients and the tests for 
equality of mean and variance of the window a single test 
statistic 


'y ;log (MAX (p ri ,10 6 ))+log P t +log p 


+ log p 

1 r l 

(9.8) 

fit of these 

tests. 

to eliminate 

seal ing 
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problems which would occur when taking the product of many 


small numbers. The p 7 's are bounded below to reduce the 

2 i 

effect of single, noisy correlation coefficients. These 

output test statistics obtained by comparing a windowed 

region of Fig. 9.3 with statistics gathered from the 

parent prototype texture sand were linearly scaled to 

produce the intensity images shown in Figs. 9.11 to 9.13. 

The results for W =W =16 are shown in Fig. 9.11. 

BIG SMALL 

Improved results for W BIG = 32 and W SMALL = are shown 

in Fig. 9.12. Results for W' = 48 and W = 16 are 

J BIG SMALL 

shown in Fig. 9.13. 

9.4 Conclusions 


In implementing the two procedures detailed in this 
chapter, large values for the test statistics Uand U 2 
were obtained even when comparing the matrices or 
correlation values measured from an extracted portion of a 
parent texture with those obtained from the entire parent 
texture. There are two basic reasons for this. First, 
the assumptions of normality are probably incorrect. 
There is little reason to believe that the distribution of 
pixels in an image is truly multi-variate normal. 
Secondly, the samples used to estimate the covariance 
matrix and the correlation coefficients are not random and 
independent. They are strongly related by their spatial 
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Figure 9.10 Input Composite 
Image 


Figure 9.11 Segmentation 
Using Indivi¬ 
dual Correla¬ 
tion Values, 



Figure 9.12 


Segmentation 
Using Indivi¬ 
dual Correla¬ 
tion Values, 



32 


Figure 9.13 


Segmentation 
Using Indivi¬ 
dual Correla¬ 
tion Values, 



48 












positions in the image and are often adjacent. The sample 
size may be adjusted to reflect this fact but no work in 
this area was done. 

In most cases, statistical tests may be considered to 
be robust - not weakened significantly - when assumptions 
used to derive them are violated. This is probably not 
true of the test for equality of two covariance matrices 
which has been called "frail at best" [47], For these 
reasons, the methods presented may be considered 
innovative but not necessarily statistically sound. 

In spite of this fact, the results show that the 
method using the test for equality of covariance matrices 
was superior to the method involving multiple tests of 
individual correlation coefficients. This could be due to 
the method used to combine the multiple tests. 

The pictorial results of this chapter indicate the 
usefulness of correlation values in texture 
identification. The methods are not intended to be 
improvements over existing segmentation identification 
techniques. An adequate number of discrimination 
techniques have been proposed by researchers already. The 
discussion is intended to apply a specific texture 
synthesis model which was based on second-order statistics 

to the texture identification problem and this was done 
successfully. 


205 








CHAPTER 10 


SUMMARY 

10.1 Introduction 

In this chapter, the results of each of the texture 
synthesis methods presented in this thesis are discussed 
and compared. This summary will clarify the methods and 
analyze the favorable and unfavorable characteristics of 
each. 

10.2 Tabulation and Discussion of All Models 

A complete synopsis of the synthesis methods 
presented in this thesis is effectively contained in 
Table 10.1. Listed in this table are the eleven methods 
of texture generation and simulation presented in 
Chapter 2 through Chapter 7. The first three columns of 
the table state the location of the text and figures 
(output) associated with each and also the key figure or 
equation number which identifies the method in a simple, 
straightforward manner. 

The next three columns are used to evaluate the 
complexity of statistics collection, statistics storage 
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and the generating process as presented in Fig. 1.1. The 
computational requirements are approximate and are 
indicated in CPU time of a DEC KL-10. Naturally, these 
numbers are relative to the processor used. The storage 
requirements are given in terms of full-word processor 
locations needed for statistics storage. In some cases, 
this storage could be reduced by packing more than one 
number (especially integers) into one 32- or 36-bit word 
but that was not done here. For example, four 8-bit 
integer values will fit into one 32- or 36-bit word. 

The seventh column provides a relative measure of the 
quality of the texture simulation on a zero to ten scale. 
A value of 5 indicates a good or reasonable simulation. 
As synthesis evaluation is a nebulous process so is the 
assignment of relative merit. The assessment of results 
is internal to this thesis as there is little synthesis 
work in the general literature. 

The last column of the table contains important 
general comments on each method. 

The one-dimensional binary generation method of 
Chapter 2 permitted study of visual response to changes in 
the probability distribution of texture. Although the 
method was not useful for natural texture simulation, it 
did lead to the models of later chapters. 
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The N-gram model of Chapter 3 was an extension of the 
Chapter 2 model applied to two-dimensional natural 
textures. The results were very good even with the severe 
constraint imposed by the upper limit on the number of 
pixels allowed in the generation kernel. With an increase 
in the complexity of the collection process this model was 
extended (see section 3.9) to allow both a larger number 
of pixels to be in the kernel and an increase in kernel 
range. The extension produced better simulations of 
structured textures. Finally, in section 3.10, the binary 
linear model, which was used to determine the contents 
(shape) of the generation kernel, was used to generate 
binary textures. The textures generated using this model 
were nearly equal in quality to those of the more complex 
and storage-consuming N-gram model. The N-gram model of 
Chapter 3 uses a generation kernel whose contents (shape) 
depends on the linear model. Therefore, the number of 
computations required in the statistics collection portion 
of the N-gram model necessarily includes computations of 
the linear model. However, in some cases, points which 
lie far from the kernel eye can be neglected in the N-gram 
model as only the best few are used due to storage 
limitations. On the other hand, such points should be 
included in the linear model therefore a larger 
neighborhood surrounding the kernel eye should be used in 
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the estimation of the linear model. 

Realizing the power of second-order statistics, a 
method of reproducing Nth order statistics using algebraic 
reconstruction was presented in Chapter 4. The approach 
proved to be academic as the number of iterations leading 
to a solution yielding adequate synthesis results was very 
large and required much storage and computation. The 
simulations were also slightly less appealing than the 
methods of Chapter 3. 

In Chapter 5, the linear autoregressive model of 
Chapter 3 was applied to 256-gray-level imagery. The 
results were good considering the vast reduction of 
information caused by the statistics collection process. 
Slightly better results were obtained by allowing the 
model to contain cross terms but the resulting complexity 
suggests that the change in texture quality is not worth 
the added effort and computational expense. Using 
non-gauss ian, non-stationary noise in the model (see 
section 5.7) produced markedly better results but with a 
requirement of slightly increased storage. 

The skip-generate method of Chapter 5 may be used to 
improve the simulation of textures having a coarse 
structure. The model produces results equal in quality to 
the linear autoregressive model of Chapter 5 while 
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requiring fewer computations in the collection process. 
The piecewise linear autoregressive model presented later 
in Chapter 6 promises an analytically superior fit but 
produced results which were not appealing enough to 
warrant additional effort. 

The field-definition model presented at the end of 
Chapter 6 is useful when generating textures composed of 
subtextures. The idea of defining fields for texture 
generation is a powerful approach for both statistics 
collection and texture synthesis but the slow response of 
the autoregressive model to boundaries produced results 
less appealing than expected. 

The best-fit model of Chapter 7 represents a 
brute-force approach to texture synthesis. Though 
computationally demanding, the final results show that 
excellent texture simulations can be generated using 
complete statistics from a relatively small neighborhood. 
The problems with a small neighborhood are seen in the 
simulation of regular textures such as raffia where the 
size of the primitives in the texture is much greater than 
the window used in the best-fit calculation. 

Combining the method of choosing a kernel shape using 
a linear model (discussed in Chapter 3) with the 
skip-generate and best-fit models would result in a very 
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powerful, but extremely computationally-demanding texture 
synthesis method. The effort required to generate 
textures by such a complex combination of methods would 
not be required in many cases where simpler models will 
produce adequate results. The complexity of the model 
depends on the texture being simulated. 

The ideas presented in Chapter 8 allow the 
introduction and removal of texture non-homogeneities. In 
Chapter 9, texture identification methods using the 
statistics employed in model parameter calculation of 
earlier chapters were proposed. Although they do not 
discriminate textures as well as the methods of other 
researchers they do illustrate the application of 
synthesis models to texture identification problems. 

10.3 Suggestions for Future Study 

Many approaches to texture analysis have been carried 
out in the frequency domain by using transform and 
filtering techniques. Texture synthesis could be carried 
out in such a domain or on frequency-filtered image. For 
example, numerous textures could be generated in 
non-overlapping frequency planes and then added together 
to obtain a final texture synthesis. However, each of 
these planes is probably interdependent and a simple 
generation with summation is probably not possible. 
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Another future approach might employ to a greater 
extent those statistics useful for texture identification 
and discrimination. In most cases, these measurements are 
not readily suitable to a synthesis process but with 
careful study, many could possibly be used in such a 
manner. Still, there is little evidence thus far to 
indicate that statistics useful for texture identification 
will be useful for texture synthesis. 

An area which deserves immediate attention involves 
preprocessing of texture by convolution. Noise filtering 
and smoothing could be useful in improving model parameter 
estimates. Also, a deconvolution processing of images 
synthesized to resemble convolved textures might produce 
excellent simulation results. 

Along this line, another synthesis method should be 
considered. With any model there is always unexplained 
variance which the model fails to account for when it is 
applied to the original parent texture. Ideally, this 
unexplained variance would be merely noise but this is 
rarely the case in practice. It is possible to develop 
additional models which could be used to explain (or 
simulate) this previously-unexplained variance. Combining 
these new models with the original model would result in 
an improved synthesis method. 
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ti .ally, texture is often produced by the shading 
effect of a light source on a three-dimensional object. 


Therefore, to generate a texture, a three-dimensional 
relief could be formed then shaded using current, 
commonly-used graphics techniques. Such an approach could 
be worthwhile in some cases but in most, synthesizing 
three-dimensional relief (which is the same as generating 
height information over a two-dimensional grid) is 
equivalent to generating intensity information over a 
two-dimensional grid. 

10.4 Conclusion 

Many natural textures are generated using a variety 
of methods presented in this thesis. The quality of the 
natural texture simulations depends on the amount of 
computation and storage used in each process. Many 
textures were adequately simulated using simple models 
thus providing a potentially great data compression for 
many applications. Others required more extensive 
computation to synthesize visually pleasing results. 
Thus, as might be expected, the success of any synthesis 
method is highly dependent on the texture itself. When 
examining the results of any method the characteristics of 
both the model and the textures used must be considered'. 
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It would be unwise to believe that all textures could 
be generated using any single approach, especially one 
which promises to compress texture information to a 
handful of numbers. Yet this is precisely what has been 
attempted in the texture synthesis work of this thesis. 

It is important to note the power and complexity of 
each synthesis method of this thesis. Many textures can 
be simulated well using simple models such as the 
autoregressive model if the model is carefully 
constructed. Improvements in texture simulation were made 
by modifying these models and allowing them to become more 
complex and use more information in the generation 
process. Other textures require more complex models such 
as the best-fit model of Chapter 7. The shortcomings of 
each method will constantly indicate where future work can 
be done. 
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APPENDIX A 


GENERATING BINARY TEXTURE PAIRS POSSESSING 
EQUAL SECOND-ORDER STATISTICS 

Texture pairs may be generated which have equal 
second-order statistics but are visually d iscriminable. 
Let G (V^, V 2 > v 3 ' V 4 ^ represent the probability of 
generating a 0 after the pixels ^j' V 2' V 3' V 4 a ^ on< 3 a l* ne 
in a one-dimensional texture (a) and G 1' V 2' V 3' V 4^ * n 
textures (b). (This is a slight change of notation from 
Chapter 2 as (a) and (b) will be texture number indices in 
this Appendix.) Define 

V- = |l“V i | (A.1) 

where V^e{ 0,1). 

The restrictions used to generate textures with equal 
second-order statistics, p a ( v i» v j) = p b^ V l ,V j^' * n Chapter 
2 may be stated as 

G a (V l' V 2' V 3'V = G a (V l' V 2 ,V 3' V 4 ) (A.2) 

WVW = <fl - 3) 
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and 


wwv = 


G b * V 1' V 2' V 3 ,V V 


(A.4) 


Equations (A.2) and (A.3) must also hold for texture 
(b). Combining Eq. (2.14) and (2.24) for a 4-gram system, 
it can be shown that 


v v 

P a (V 


1 


1 



' V 3'V 


' V 3'V 


P 

a 


(V 1' V 2 ,V 3' V 4 ) 
(V 1' V 2' V 3' V 4 ) * 


(A.5) 

(A.6) 


Combining the above restrictions yields two further rules, 

V V' 

[v 5 +(-l) 5 G a (v 1 ,v 2 ,v 3 ,v 4 ))=[v£+ (-1) 5 G b (v 1 ,v 2 ,v 3 ,v 4 )] (A - 7) 

V v' 

[V 5 +(-l) 5 G a (V 1 ,V 2 ,V 3 ,V 4 ))=[V'+(-l) 5 G a (V 1 ,V 2 ,V 3 ,V 4 )] 

v' 

= [V'+(-l) 5 G a (V ; [,V 2 ,V 3 ,V 4 )] ( f* 8) 
These restrictions yield textures having the equalities 


p a (° )= P b (0) = P a (l) = P b (l) • (A.9) 

As no closed-form solutions to Eq. (2.14) and Eq. (2.24) 
are easily obtained the proofs verifying Eq. (A.5) and 
(A.9) are very complex. However, these properties may be 
illustrated by generating numerous examples where 
Eq. (A.2), Eq. (A.3) and (A.4) hold. 
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Assuming the above conditions, a proof then follows 
to show that 
restrictions hold: 


P a (Vi,Vj) = Pfc)( v l'Vj) if the above 


P a (V l' V j> ' Z P a<V V 2.Vl'V 

V. 

1 


i^i / j 


-ZvwwJ 

=Z P b ( W V 3' V 4> « 'V | - 1) ' 5s b (V W' V k-3' V t-2’ , k-l 11 


" ’V'-D \ l, H' , k-J'V2'Vl" 

k=5 

j 
n 

k=5 


=y^P h (v 1 ,v^,v 3 ,v 4 ) . 


(A.10) 


M 


v: 


n '{lv;^+(-D 3k+2 G^(v,,. ,)] 


k=5 


3k-2 


b 3k-2 3k-l 3k 3k-l 


tV 3k-3 +(_11 3G b (V 3k-l ,V 3k' V 3k-l’ V 3k-2 )1 


|V 3k + 4 + <- X 1 3k+4<3 b <V 3k' V 3k-1 ' V ik-2' V 3k + 3 
v: 


»} 


^ [V j+l-H +< - 1 > j + 1 A < '' j - t -3' V 5-*-2' V j-*-l' V j-t n 
where r = M0D(j-l,3). It is very important to note that 
the variables in the product expression match 
successively. That is, V^k-l anc ^ V 3k+2 a ^ t |ave the prime 
notation in each sub-expression. At this point, there are 
three possible cases: 




0 


Case 1, r = 


In this case, the proof is completed by a series of 

change of variables as summing over is the same as 

summing over V'. The remainder of the proof becomes 

V, 


E 3 

P fo( Vi » V 2 » V 3 ' V 4 )* n [V k +(-l) G b (V k-4 ,V k-3 ,V k-2 ,V k-l 
k=l , k 


)] 

(A.11) 


- vw 


Case 2, r = 1 


In this case, the remainder of the proof becomes 

j- 1 v k 

P b (V l' V 2' V 3' V 4 } n ^ v k + <~ 1 ) G b (V k-4' V k-3' V k-2' V k-l )1 * 

k—1 


V. 


.[V' + (-n 3 G b (V j _ 4 ,V'. 3 , Vj _ 2 .V j . 1 


)] 


j-1 


(A.12) 


P b (V l' V 2' V 3' V 4 ) lr n 1 IV k +( “ 1) G b (V k-4' V k-3' V k-2'Vl )] 

k=l 

.[vr + (-l) JG b (v j . 4 ,v.. 3 ,v j . 2 ,v.. 1 )l . 


Thus, in this case P a (V lr Vj) = P b (V lf Vj). This implies 
that p a (0 1 f0j) = p b (0 1 flj); p a (0j/ lj 1 = p b (0 l'°j ); 

= P b (l r lj)» p a (li»lj) = P b a v 0].Vie know P a (O r lj) - Wj) 
and P b (O r lj) = P^l^Oj) and P (0) = P (1) = 0.5. So P a (V lf V.j) 
= P b (Vi,Vj) = 0.25 for all and Vj . 
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Case 3, r_ = 2_ 

In this case, as in the case where r = 0, 
remainder of the proof is straightforward, requiring only 
a change of variables. 
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APPENDIX B 

GOODNESS-OF-FIT TEST 
FOR N-GRAMS 

We may think of the one-dimensional binary texture 
generation process as an experiment which should yield 
certain N-grams given certain generation parameters. 
These N-grams or Nth-order densities may be determined 
analytically by combining Eq. (2.14) and Eq. (2.24) of 
Chapter 2. However, as with any Monte-Carlo simulation, 
the analytic parameters may not agree exactly with the 
statistics or estimates of those parameters based on 
observations from the simulation or experiment. 
Statistical tests may be used to determine whether the 
statistics match the parameters to the extent expected 
given a given random sample. 

The most common statistical test for this purpose is 
the "goodness-of-fit" test involving the chi-square 
distribution. This test is used when we want to compare 
an observed distribution with the corresponding values of 
a theoretical distribution. The test is based on the 
statistic 
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2 N .2 

1=1 1 


(B.l) 


The and the e^ are the observed and expected 

frequencies of a certain N-gram pattern in our experiment. 

The sampling distribution of this statistic is 

N-l 

approximately the chi-square distribution with 2 
degrees of freedom. (The constraints on the system of 
N-grams yield this number of degrees of freedom.) A 
chi-square distribution with n degrees of freedom is given 
by 


2 (V 2 )^2 (n-2) _y 2 /2 

fix ) - - JL - ! - l X ' (B.2) 

2 n/2 r(|) 

The approximation is close provided e >5. 

As a first step in the testing process, the complete 
set of N-grams, P ^ V l ' V 2 ' * • • , v n ) i- 5 computed using 
Eq. (2.14) and Eq. (2.24). Their corresponding 

statistics, P(V^,V 2 ,...,V N ) are computed from the 
generated texture by counting the number of occurrences of 

each pattern (V^,...,V N ) and dividing by the total sample 

/\ 

size, M. Then, f i = P(V^,...,V N )*M and 

e^ = P (V^, .. . ,Vjj) *M. For large N, the number of degrees 
of freedom is greater than 30 and a normal distribution 
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is 


qui te 


accurate. 


The 


expression 


approximation 


- 


/2N-1 is approximately normally distributed as the 


standard normal distribution. The tables listed in [18] 
may also be used for large N. Additional details 


concerning the goodness-of-fit test may be found in 
[45,63-65]. 


For each of the textures in Chapter 2, this test was 
applied to each texture half and the hypothesis that the 
expected and measured N-grams are statistically equal was 
accepted at a a= 0.01 confidence level for N = 1,...,10. 
Actually, these tests are not independent in a statistical 
sense as the information content for each N overlaps. But 
this poses no violation of assumptions of the statistical 
test. However, all possible patterns in the generated 
sample may be a violation of the independence of samples 
assumption. Still, a low chi-square value indicates a 
good fit of the data to predicted parameter. 
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APPENDIX C 

PSEUDO-RANDOM VARIATE GENERATION 

Two algorithms were used in this thesis to generate 
uniformly-distributed pseudo-random numbers. Both are 
modifications of the multiplicative congruence method 
[48,66]. These uniformly-distributed variates may be 
transformed to normally-distributed deviates using the 
inverse normal probability distribution function. 

The first algorithm is the Lehner Pseudo-Random 
Number Generator [67]. The general form of this generator 
is - 

X n+1 = KX n (M0D 231_1) (C.l) 

where 

K = 14 29 (MOD 2 31 -1) = 630360016 10 

The resulting number is basically a 31-bit pattern which 
is in the form of a floating-point number in the range 
(0,1). The algorithm can be written to avoid division by 
2 31 -1 [68 1 . 

The second algorithm used to generate uniformly 
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distributed variates is 

X = 7 5 X (MOD 2 31 -1) . 

n+1 n 


(C. 2) 


i t 
in 

to 
is 

x = -i- fV t2 / 2 at . (c - 3) 

✓ 5 ? J 

—oo 

The algorithm for accomplishing this is described in Ref. 
[57] and is briefly outlined here. 

The basic interval (0,1) is divided into 4 segments. 
In each segment the inverse of the Gaussian integral, 
invgauss(X), or an integral of a similar form is 
approximated by a minimax rational function [55,72], The 
approximations for the final 3 segments (comprising the 
interval Xe{(0,0.075) (0.925,1)}) are functions of the 
transformed variable 

W = V-log e(1-X 2 ) . (C.4) 

This transformation of the variable improves the 


The resulting integer is multiplied by 2 to convert 
to a floating-point number. This generator is reported 
Refs. [69-71]. 

These uniform deviates may be transformed 
normally-distributed deviates if desired. this 
accomplished by computing the inverse of the integral 
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efficiency and stability of the approximation [SO], The 
rational functions in these intervals are of the form 

( C^W+C-W +C 3 W 

c n + -T- 3 

d Q +d 1 W+d 2 W +W 

For the remainder of the (0,1) interval, 
xe{(0,075, 0.925) }, the function is of the form 

invgauss(X) = 

(C . 6 ) 

(Z+Z-(b 0 +a 1 .Z 2 /(b 1 +Z 2 +a 2 /(b 2 +Z 2 +a 3 /(b 3 +Z 2 ))))) *SGN(Z) 

where z = 11 —2 x |. The constants a^, bj_, c ^ and d^ have 

specific values found by.solving the minimax problem and 
are given in Ref. [57]. 


)■ 


(C. 5) 
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